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1.  INTRODUCTION  AND  OVERVIEW 


The  Tropospheric  Electromagnetic  Parabolic  Equation  Routine  (TEMPER)  has 
proven  to  be  a  useful  computational  tool  for  modeling  propagation  through  the 
troposphere.  The  routine  uses  a  scalar  paraxial  approximation  to  the  vector 
electromagnetic  wave  equation.  The  resulting  equation  can  be  solved  by  a 
marching  procedure;  given  an  initial  field  disU'ibution  and  a  description  of  the 
medium,  the  vertical  field  distribution  can  be  calculated  at  arbitrary  range.  Standard 
atmospheric  models  or  experimentally  derived  data  are  used  to  quantify  the 
macroscopic  features  of  the  index  of  refraction.  When  sufficient  atmospheric  data 
are  available,  TEMPER  yields  generally  good  agrees,  ent  with  field  experiments. 
(Dockery*  and  Kuttler  and  Dockery ,2  and  the  references  therein,  provide  a  more 
detailed  discussion  of  this  method.) 

In  the  actual  atmosphere,  small-scale  fluctuations  in  the  index  of  refraction  are 
superimposed  upon  the  large-scale  features.  These  fluctuations  are  random  and 
highly  chaotic  and  are  caused  by  atmospheric  turbulence.  Since  the  physical  scale 
of  the  turbulence  can  range  from  perhaps  100  m  to  as  little  as  1  mm,  it  is  clearly 
impractical  to  conduct  environmental  measurements  with  sufficient  resolution  to 
characterize  completely  a  medium  that  might  extend  for  many  kilometers. 
Although  the  magnitude  of  these  perturbations  may  be  small,  they  can  serve  to 
focus  and  defocus  the  propagating  field.  At  sufficiently  long  ranges,  the 
cumulative  effect  can  be  significant.  Consequently,  it  is  desirable  to  include  the 
statistical  nature  of  these  random  perturbations  in  the  propagation  model. 

This  report  summarizes  the  progress  to  date  of  an  ongoing  effort  to  incorporate 
random  refractivity  fluctuations  in  TEMPER.  In  the  proposed  approach, 
appropriate  spectral  models  for  the  surface  layer  turbulence  are  developed. 
Individual  realizations  consistent  with  the  desired  spectrum  are  generated 
numerically.  The  field  is  then  propagated  through  the  medium  using  the  parabolic 
equation/split-step  algorithm.  By  averaging  across  an  ensemble  of  realizations,  the 
statistical  properties  of  the  propagated  field  are  estimated.  The  numerical 
simulations  presented  in  this  report  display  excellent  agreement  with  existing 
theoretical  expressions  for  canonical  wave  propagation  problems.  These  solutions 
provide  a  benchmark  for  evaluating  the  numerical  algorithm. 

Although  the  parabolic  equation  method  has  been  widely  used  in  optics  and 
acoustics  to  study  wave  propagation  through  random  media,  electromagnetic 
propagation  through  the  troposphere  presents  a  unique  set  of  problems.  Unlike  the 
situations  typically  encountered  in  optics  and  acoustics,  electromagnetic 
propagation  is  usually  governed  by  weak  scattering  theory.  This  difference  in 
theoretir  d  approach  affects  the  statistical  quantities  that  are  estimated  in  the 
numcrica  simulations.  Unlike  the  specha  in  analogous  ocean  acoustics  problems, 
multidimensional  special  representations  for  atmospheric  turbulence  are  not 
separable  into  the  products  of  simple,  one-dimensional  spectra.  In  addition,  unlike 


’G.  D.  Dockery,  Description  and  Validation  of  the  Electromagnetic  Parabolic 
Equation  Propagation  Model  (EMPE),  JHU/APL  FS-87-152,  Applied  Physics 
Laboratory,  Laurel,  Md.  (1987). 

R.  Kuttler  and  G.  D.  Dotkery,  "Theoretical  Description  of  the  Parabolic 
Approximation/Fourier  Split-Step  Method  of  Representing  Electromagnetic 
Propagation  in  the  Atmosphere,"  Radio  Sri  26,  381-394  (1991) 
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the  homogeneous  turbulence  of  the  free  atmosphere  usually  encountered  in  optics, 
surface  layer  turbulence  at  low  altitudes  is  inhomogeneous  and  possibly 
anisotropic.  These  features  and  others  have  a  significant  effect  on  the  way  in 
which  turbulence  is  modeled  and  coupled  into  a  parabolic  equation  routine. 

Section  2  briefly  reviews  the  properties  of  homogeneous  isotropic  turbulence 
and  the  associated  Kolmogorov  spectrum.  Existing  statistical  models  for 
atmospheric  turbulence  are  described  and  modified  spectra  useful  in  wave 
propagation  studies  are  considered.  The  effects  of  the  turbulent  inner  and  outer 
scales  on  the  spectrum  are  noted.  The  structure  constant  of  tu  rbulence  is  related  to 
the  macroscopic  mean  refractivity  profile.  The  height-dependent  nature  of  the 
turbulence  for  low  altitudes  is  modeled,  and  the  profile  used  by  Levy  and  Craig^  in 
their  parabolic  equation  simulations  is  presented  as  a  numerical  example.  The 
advantages  and  limitations  of  spectral  modeling  for  simulating  turbulence  are 
discussed. 

Section  3  begins  by  reviewing  the  classical  theories  for  wave  propagation 
through  random  media.  Both  weak  and  strong  scattering  theories  are  discussed. 
Whereas  strong  scattering  theory  is  concerned  with  the  statistical  moments  of  the 
total  field,  weak  scattering  theory  typically  deals  with  the  moments  of  either  the 
log-amplitude  or  phase  of  the  field.  Existing  results  for  plane  wave  and  point 
source  propagation  are  summarized.  The  validity  of  using  two-dimensional 
simulations  such  as  TEMPER  to  model  three-dimensional  propagation  through  a 
random  medium  is  studied  for  a  plane  wave.  It  is  shown  that  although  two- 
dimensional  modeling  often  suffices  for  calculating  the  moments  of  the  field,  it 
may  be  inadequate  for  calculating  the  log-amplitude  fluctuations.  A  numerical 
example  is  given,  illustrating  the  differences  between  two-  and  three-dimensional 
propagation. 

Section  4  details  how  random  fluctuations  can  be  included  in  TEMPER.  The 
Unnsverse  spectrum  is  shown  to  be  the  necessary  statistical  quantity  for  describing 
the  medium.  Analytical  expressions  for  the  tfansverse  specmim  of  turbulence  are 
derived  and  presented  in  terms  of  special  mathematical  functions.  Simplified 
expressions  valid  at  microwave  frequencies  are  then  deduced  and  compared  to  tliose 
of  Levy  and  Craig.^  Numerical  realizations  of  the  transverse  process  are  generated 
and  the  multiscale  properties  of  turbulence  are  observed.  Results  of  preliminaj  y 
wave  propagation  studies  are  presented.  The  three-dimensional  effects  in  the  log- 
amplitude  fluctuations  predicted  in  section  3  are  confirmed  by  numerical  experi¬ 
ments.  Good  agreement  with  three-dimensional  theory  is  obtained  for  the  second 
moment  of  the  scattered  field. 

Finally,  section  5  summarizes  the  main  results  of  this  study.  A  program  for 
upgrading  TEMPER  by  including  random  refractivity  fluctuations  is  proposed. 


^M.  F.  Levy  and  K.  H.  Craig,  "Millimetre- Wave  Propagation  in  the  Evaporation 
Duct,"  in  Atmospheric  Propagation  in  the  UV,  Visible,  IR  and  MM-Wave  Region 
and  Related  Systems  Aspects,  AGARD  Conf.  Proc  454,  Neuilly-sur-Seine,  France, 
pp.  26.1-26.10  (198'9). 
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2.  STATISTICAL  MODELS  FOR  ATMOSPHERIC 
TURBULENCE 


In  contrast  to  steady  laminar  flow,  turbulence  is  chaotic  and  characterized  by 
random,  highly  rotational  motion.  Atmospheric  parameters  such  as  wind  speed, 
temperature,  humidity,  and  pressure  can  all  exhibit  turbulent  fluctuations.  Since 
some  of  these  quantities  determine  the  index  of  refraction,  an  electromagnetic  wave 
propagating  through  a  turbulent  medium  will  also  acquire  random  variations.  To 
predict  the  statistical  nature  of  the  field  fluctuations,  it  is  necessary  to  describe  the 
statistical  properties  of  the  medium  through  which  it  propagates. 

This  section  will  emphasize  the  use  of  existing  spectral  models  for  turbulent 
fluctuations  in  the  index  of  refraction.  The  spectral  modeling  approach  is  attractive 
for  several  reasons.  For  the  theorist,  existing  wave  propagation  theories  are 
usually  formulated  in  terms  of  the  spectrum  of  the  random  medium.  For  the 
experimentalist,  spectra  can  often  be  quantified  by  a  limited  number  of  observable 
parameters.  Spectral  models  are  also  easy  to  implement  numerically;  to  simulate 
realizations  of  an  assumed  Gaussian  random  process,  all  that  is  required  is  a  random 
number  generator  and  a  fast  Fourier  transform  (FFT)  routine.  The  approach  can  be 
adapted  to  temporal  spectra  and  to  spatial  spectra  of  any  number  of  dimensions. 
Finally,  of  particular  relevance  to  upgrading  TEMPER,  the  spectral  approach  has 
already  been  extensively  applied  in  conjunction  with  the  split-step  algorithm  to 
study  acoustic  propagation  through  a  random  ocean. 

2.1  LOCALLY  HOMOGENEOUS  ISOTROPIC  TURBULENCE 

Turbulence  is  a  cascade  process  whereby  large-scale  eddies  are  broken  into 
successively  smaller  sizes  until  they  are  eventually  dissipated  in  the  form  of  heat. 
The  outer  scale,  Lq,  for  turbulence  in  the  free  aunosphere  typically  ranges  from  10 
to  100  m.  The  inner  scale,  t  ^Jso  called  the  Kolmogorov  microscale,  is  on  the 
order  of  millimeters  or  centimeters.  Eddies  larger  than  Lq  feed  energy  into  the 
turbulence.  Between  the  two  extremes,  the  eddies  are  said  to  be  in  the  inertial 
subrsrge,  where  the  kinetic  energy  of  the  eddies  dominates  the  dissipative  effects. 
Structures  less  than  are  quickly  dissipated  owing  to  viscosity. 

Fully  developed  turbulence  is  represented  by  the  Kolmogorov  spectrum.  Since 
fully  developed  turbulence  is  isotropic,  the  three-dimensional  spectrum  of  the  index 
of  refraction  can  be  written  as  a  function  of  a  single  variable  K,  where 

=kl+ky  +kj  .:nd  kj^.ky,  andjfej  represent  the  spatial  frequency  components 
in  the  different  directions.  As  shown  in  Figure  2.1,  the  spectrum  can  be  divided 
into  three  regimes.  No  general  expression  applies  for  spatial  frequencies 
0<K  <2nl  Lq,  where  energy  is  input  into  the  turbulence.  Within  the  inertial 
subrange  2nl  Lq<K  <2nUQ,t)\c  spectrum  is  described  by  a  power  law  and  goes 
like  In  the  dissipation  range,  the  spectrum  falls  rapidly  to  zero. 

In  practice,  meteorological  parameters  are  often  measured  as  a  function  of  a 
single  spatial  variable  from  which  one-dimensional  spectra  V{K)  are  estimated. 
Since  for  an  isotropic  process  the  three-dimensional  spectrum  is  also  a  function  of 
a  single  variable,  the  two  spectra  for  the  refractive  index  n  are  related  by  (see  Ref. 
4,  p.  517) 


‘*A.  Ishimani,  Wave  Propagation  and  Scattering  in  Random  Media,  Ac.idemic  Press, 
New  York  (1978). 
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Figure  2.1  Kolmogorov  spectrum  for  turbulence.  Spatial  frequencies  0  <  K  < 
2n/Lo.  where  Lq  Is  the  outer  scale,  constitute  the  energy  input  region.  No  general 
form  for  the  spectrum  exists  in  this  range.  The  inertial  subrange  is  bounded  by 
2nlL()<K <2nUQ,  where  Iq  is  the  inner  scale  of  the  turbulence.  Within  this 
region,  the  three-dimensional  spectrum  is  proportional  to  The  spectrum 

rapidly  decreases  in  the  dissipation  region,  where  K  >  27iUq. 
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(2.1) 


Sn{K)  = 


2nK  oK 


Consequently,  V„  is  proportional  to  in  the  inertial  subrange.  The  transverse 
spectrum,  which  is  also  one-dimensional  but  obeys  a  different  power  law,  is 
defined  in  section  4.1.  It  will  be  shown  that  the  transverse  spectrum  and  not  VniK) 
is  of  relevance  for  including  turbulent  fluctuations  in  TEMPER  simulations. 

Usually,  a  spectrum  is  related  to  a  correlation  function  by  a  Fourier  transform 
pair  through  the  Wiener-Khinchin  theorem.  For  this  transform  pair  to  exist,  the 
process  must  be  statistically  homogeneous.  Statistical  homogeneity  is  an 
extension  of  the  familiar  time  series  concept  of  wide  sense  stationarity  to  multiple 
spatial  dimensions.^  A  statistically  homogeneous  process  will  have  a  constant 
mean  and  its  autocorrelation  will  depend  only  on  the  distance  between  the  points 
being  correlated  and  not  on  the  absolute  position.  Unfortunately,  atmospheric 
turbulence  is  not  statistically  homogeneous,  since  both  of  the  mean  characteristics 
will  fluctuate  in  space  and  the  correlation  function  will  depend  on  position.  The 
lack  of  homogeneity  is  reflected  in  the  Fourier  domain  by  the  spectrum  being 
undefined  for  frequencies  0  <  K  <  IuILq.  To  circumvent  these  difficulties, 
Kolmogorov  extended  the  concept  of  stationarity  to  define  a  process  that  is  locally 
homogeneous.  He  noticed  that  although  the  wind  velocity  is  not  homogeneous, 
the  difference  between  the  wind  velocities  taken  at  points  less  than  Lq  apart  will 
behave  like  a  homogeneous  process.  This  locally  homogeneous  process  is  called  a 
structure  function.  Denoting  the  three-dimensional  position  vector  by  r,  the 
structure  function  of  the  index  of  refraction  is 


0„(r)  =  ^|/i(r, +r)-«(ri)f^. 


(2.2) 


where  the  angle  brackets  indicate  an  ensemble  average.  If  the  process  is  also 
isotropic,  then  D„  =  D„(r),  where  r  =  Irl.  For  an  isotropic,  locally  homogeneous 
Kolmogorov  process,  the  structure  function  is  given  by 


0„(r)  =  C2r^ 


for  LQ»r»  to 


(2.3a) 


D„ir)  =  C^t^{r/to)^  for  r«to. 


(2.3b) 


where  C„  is  the  structure  constant.  Equation  2.3a  is  a  statement  of  the  well-known 
"Kolmogorov  two-thirds  law"  that  permeates  the  turbulence  literature.  The 
structure  function  is  not  explicitly  defined  for  distances  greater  than  Lq,  which  acts 
like  the  correlation  length  of  the  medium.  For  distances  greater  than  the 
correlation  length,  however,  /i(ri-i-r)  and /i(ri)  in  Equation  2.2  will  become 


^A.  J.  Devaney,  "The  Inverse  Problem  for  Random  Sources,"  J  Math  Phys.  20, 
1687-1691  (1979). 
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essentially  uncorrelated  and  the  structure  function  will  become  related  to  the 
variance  of  the  process.  Finally,  the  structure  function  of  an  isotropic  process  is 
related  to  the  spectrum  by 


Dn(r) 


=  8rrr[l-( 

Jo 


(Kry^sin(Kr)]S„(K)K^dK 


(2.4) 


(see  Ref.  4,  p.  524).  Note  from  Equation  2.4  that  spatial  frequencies  K  <  InlLo 
make  minimal  contributions  to  the  integration  for  r«Zoand  reasonable  spectra; 
consequently,  the  structure  function  and  the  spectrum  can  be  related  by  a  transform 
pair  even  when  is  not  well  known  for  low  frequencies. 

2.2  MODIFIED  SPECTRAL  REPRESENTATIONS 

Figure  2.1  shows  three  separate  spectral  regimes  for  turbulence  with  the  energy 
input  regime  undefined.  For  mathematical  convenience,  it  is  often  useful  in  wave 
propagation  theory  to  assign  a  spectral  shape  to  the  input  regime  and  then  to 
combine  all  three  regimes  into  a  single  function.  The  resulting  expression  is  the 
so-called  von  Karman  spectrum  (see  Ref.  4,  p.  337,  and  Ref.  6,  p.  5), 


S„(K)  =  Q.Q3K^(k^  +  / Ki),  (2.5) 


where  K„  =  5.92l Iq.  Equation  2.5  must  be  considered  an  approximate 
representation  of  turbulence,  since  the  spectral  shape  for  low  frequencies  is 
unknown.  The  von  Karman  specuum  is  a  well  behaved  function  that  has  a 
conventional  inverse  Fourier  transform  and  hence  a  function  that  can  be  interpreted 
as  the  autocorrelation:  the  von  Karman  spectrum  consequently  can  describe  a 
process  that  is  bcth  locally  homogeneous  and  statistically  homogeneous.  These 
distinctions  between  locally  homogeneous  and  stati.siically  homogeneous  processes 
are  subtle  and  are  discussed  in  further  detail  by  Ishimaru'*  and  Tatarskii.^  The  main 
point  is  that  if  the  von  Karman  spectrum  is  used  in  any  wave  propagation  analysis, 
the  resulting  solutions  for  the  statistics  of  the  field  are  at  best  reliable  only  for 
correlation  distances  less  than  the  outer  scale  of  the  medium.. 

Kolmogorov's  original  insights  were  based  on  dimensional  analysis  and  the 
assumption  that  the  turbulence  is  isotropic.  For  well  developed  turbulence  far 
from  any  boundaries,  the  isotropic  assumption  is  usually  valid.  At  low  altitudes, 
however,  turbulence  becomes  anisotropic.  For  example,  as  mentioned  in  section 
1,  a  typical  value  for  the  outer  scale  in  free  turbulence  is  100  m.  At  altitudes  less 
than  100  m,  the  horizontal  outer  scale  might  remain  unchanged,  but  the  vertical 
correlatio'i  scale  must  be  reduced.  The  Kolmogorov  spectrum  can  be  modified  to 
model  an  anisotropic  process  yielding  (Ref.  4,  p.  339) 


®B.  J.  Uscinski,  The  Elements  of  Wave  Propagation  in  Random  Media,  McGraw-. 
Hill.  New  Vork  (1978). 

’V.  I.  Tatarskii,  The  Effects  of  the  Turbulent  Atmosphere  on  Wave  Propagation, 
Keier  Press,  Jerusalem  (1971), 
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Sn(K,ky,k,)=om3c^(kltl+k^tl+k’^i]) 


-11/6 


(2.6) 


where  tx,  ty,  and  are  the  outer  scales  in  the  three  directions.  The  lowest  1(X) 
m  of  the  atmospheric  boundary  layer,  where  turbulence  is  often  anisotropic,  is 
called  the  surface  layer.*-®  Equation  2.6  is  valid  in  the  inertial  subrange. 
Presumably,  some  hybrid  of  Equations  2.5  and  2.6  could  be  developed  to  yield  a 
complete  spectral  representation  of  anisotropic  surface  layer  turbulence. 

For  simplicity,  the  isotropic  model  of  Equation  2.5  will  be  used  in  the  majority 
of  this  report.  This  is  done  with  the  understanding  that  a  parallel  development  is 
possible  for  the  anisotropic  model  should  it  prove  necessary. 

A  further  complicating  factor  is  that  the  "structure  constant"  C„  is,  in  fact,  not 
constant  for  altitudes  in  the  surface  layer.  This  is  discussed  in  further  detail  in 
section  2.3. 

2.3  STRUCTURE  CONSTANT  AND  SURFACE 
LAYER  TURBULENCE 

The  structure  constant  C„  determines  the  strength  of  the  turbulence.  In  the  free 
atmosphere,  it  ranges  typically  between  10”^  m“*^  for  strong  turbulence  and  10"® 
m"*/*  for  weak  fluctuations.  If  the  turbulence  is  approximated  by  a  homogeneous 
process  by  assuming  the  von  Karman  spectrum,  C„  can  be  related  to  the  variance 
of  the  random  part  of  the  index  refraction  fluctuations  <nj  >  and  the  outer  scale 
(Ref.  4,  p.  543): 


cl^\.9\{n})Lf\ 


(2.7) 


Clearly,  the  numerical  determination  of  the  structure  constant  for  a  given  wave 
propagation  problem  is  of  crucial  importance;  as  will  be  shown  in  a  later  section, 
the  solution  for  the  electromagnetic  field  fluctuations  will  be  specified  in  terms  of 
C^.  One  approach  to  calculating  cl  begins  by  specifying  the  index  of  refraction 
in  terms  of  the  measurable  meteorological  parameters.  A  standard  expression  for 
the  microwave  index  of  refraction  is 


„_1=  7M(p  +  48iOe/7)xlO'^,  (2.8) 


where  T  is  the  temperature  in  degrees  Kelvin,  P  is  the  pressure  in  millibars,  and  e 
is  the  water  vapor  pressure  in  millibars.  Rather  than  use  the  index  of  refraction, 
meteorologists  prefer  to  use  the  potential  index  of  refraction  and  write  an  analogous 
expression  in  terms  of  the  potential  temperature  and  the  potential  vapor  pressure 


*H.  A.  Panofsky  and  J.  A.  Dutton,  Atmospheric  Turbulence  Models  and  Methods 
for  Engineering  Applications,,  John  Wiley  and  Sons,  New  York  (1984). 

^Z.  Sorbjan,  Structure  of  the  Atmospheric  Boundary  Layer,  Prentice  Hall,  Engle¬ 
wood  Cliffs,  N.  J.  (1989). 
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(specific  humidity).  These  potential  functions  are  called  conservative  passive 
additives.  Unlike  T  and  e,  a  volume  element  of  a  conservative  passive  additive 
preserves  its  value  as  it  moves  through  space,  and  it  does  not  exchange  energy  with 
the  turbulence.  If  we  are  only  interested  in  refractivity  at  a  constant  height,  no 
distinction  is  necessary.^ ‘ 

The  index  of  refraction  is  written  as  the  sum  of  a  mean  deterministic  part  <n> 
that  may  vary  slowly  in  height  and  range  owing  to  atmospheric  stratification,  and 
by  a  small,  randomly  fluctuating  part  rif  due  to  turbulence.  By  expanding  Equa¬ 
tion  2.8,  rif  can  be  written  in  terms  of  small  perturbations  in  temperature, 
pressure,  and  humidity.  The  variance  of  the  randomly  fluctuating  part  of  the  index 
of  refiraction  can  thus  be  expressed  as  a  function  of  the  variances  of  the  individual 
meteorological  parameters  along  with  certain  cross  terms  that  include  coupling 
between  the  atmospheric  variables;  for  example,  temperature  and  humidity 
fluctuations  are  usually  correlated.  The  result  is  an  expression  for  in  terms  of 
the  presumably  measurable  structure  constants  of  the  meteorological  parameters; 
various  expressions  are  available  in  the  literature.*-**'*^  The  relative  importance 
of  the  various  terms  depends  on  the  atmospheric  conditions,  the  altitude  of  interest, 
and  the  frequency  of  the  probing  wave.  The  structure  constant  at  optical 
frequencies  is  primarily  a  function  of  temperature  fluctuations,  whereas  in  the 
microwave  range  both  the  humidity  and  the  coupling  between  temperature  and 
humidity  must  be  considered. 

If  the  mean  profile  of  the  index  of  refraction  is  known  (presumably  by 
measuring  the  atmospheric  parameters),  can  be  written  in  terms  of  the  outer 
scale  of  the  turbulence  and  the  gradient  of  the  mean  profile.  Tatarskii’  gives 


cl  =  a^a'4^^(d  <n>  Idzf ,  (2.9) 


where  r  is  the  height  and  a'  is  a  correction  term.  The  constant  is  taken  to  be  a 
"universal  constant"  but  its  precise  value  is  not  well  known.  Various  estimates 
place  it  between  1.5  and  3.5.  For  the  atmospheric  boundary  layer,  Sorbjan’  uses 
1.6.  Ishimaru^  follows  Monin  and  Yaglom*^  and  uses  2.8.  In  a  series  of  papers, 
Gossard**'*^’*®  uses  3.2.  Panofsky  and  Dutton*  make  no  pretense  of  knowing 
this  "universal  constant"  to  any  more  than  one  significant  figure  and  use  3.  The 
correction  term  a'  suggested  by  Tatarskii  is  invariably  taken  to  be  1. 


’°D.  E.  Kerr  (ed.).  Propagation  of  Short  Radio  Waves,  Dover,  New  York  (1951). 

"E.  E.  Gossard,  "Clear  Weather  Meteorological  Effects  on  Propagation  at  Fre¬ 
quencies  above  1  GHz,"  Radio  Sci.  16,  589-608  (1981). 

*^E.  E.  Gossard,  "Refractive  Index  Variation  and  Its  Height  Distribution  in  Different 
Air  Masses,"  Radio  Sci.  12,  89-105  (1977). 

**S.  D.  Burk,  'Temperature  and  Humidity  Effects  on  Refractive  Index  Fluctuations  in 
Upper  Regions  of  the  Convective  Boundary  Layer,"  J.  Appl.  Meteorol.  20,  717- 
721  (1981).  2 

*‘*E.  L.  Andreas,  "Estimating  Over  Snow  and  Sea  Ici.  fron  Meteorological  Data," 
J.  Opt.  Soc  Am.  A.  Opt  Image  Sci.  5,  481-495  (198i.'’' 

’^A.  S.  Monin  and  A.  M.  Yaglom,  Statistical  Fluid  Mechanics,  MIT  Press,  Cam¬ 
bridge,  Mass.  (1971). 

*®E.  E.  Gossard.  "The  Height  Distnbution  of  Refractive  Index  Structure  Parameter  m 
an  Atmosphere  Bemg  Modified  by  Spatial  Transition  at  Its  Lower  Boundary," 
Radio  Sci.  13,  489-500  (1978). 
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Even  if  the  mean  profile  is  well  known  and  a  reasonable  value  for  c?-  is  selected, 
it  is  still  difficult  to  calculate  €„  by  means  of  Equation  2.9,  owing  to  its 
dependence  on  the  outer  scale  Lq.  As  noted  by  Van  Zandt  et  al.,*^  the  outer  scale 
is  notoriously  difficult  to  measure,  although  the  outer  scale  for  turbulence  in  the 
free  atmosphere  is  usually  taken  to  be  between  10  and  100  m.  By  examining  the 
backscatter  from  free  atmosphere  turbulence.  Van  Zandt  estimated  that,  on  average, 
Iq  =  10  m. 

For  altitudes  less  than  about  100  m,  the  turbulence  can  no  longer  be  assumed  to 
be  locally  homogeneous  and  isotropic.  A  first-order  expression  for  the  outer  scale 
as  a  function  of  height  z  is 


(2.10) 


where  is  von  Karman's  constant,  0  is  the  normalized  wind  sheer,  and  (p£  is  the 
normalized  dissipation.  The  von  Karman  constant  is  traditionally  taken  to  be  0.4, 
with  values  measured  in  the  wind  tunnel  and  the  atmosphere  ranging  from  0.35  to 
0.43  (see  Ref.  8,  p.  122,  and  Ref.  9,  p.  73).  For  a  neutral  atmosphere,  the 
normalized  parameters  are  equal  to  1.  Panofsky  and  Dutton*  give  formulas  for  the 
parameters  in  both  stable  and  unstable  conditions.  Combining  Equations  2.9  and 
2.10  and  setting  a'=  1  yields  the  structure  constant  in  the  boundary  layer: 


<n>  Idzf .  (2.1 1) 

2.4  NUMERICAL  EXAMPLE:  PROFILE  OF  LEW 
AND  CRAIG 

The  ultimate  goal  of  this  study  is  to  generate  numerical  realizations  of 
atmospheric  turbulence  and  then  do  wave  propagation  studies  using  the  parabolic 
equation.  A  similar  study  was  attempted  by  Levy  and  Craig.’"  In  order  to 
eventually  compare  numerical  results,  we  now  calculate  the  structure  constant  for 
the  index  of  refraction  profile  used  m  their  study. 

Referencing  a  report  by  Battaglia,'*  and  after  converting  from  modified 
refractivity  to  the  index  of  refraction.  Levy  and  Craig  use  as  their  profile 


<n(z)>-l  =  10‘^(a-10^/a^)  z- ad  10“^  log (-^j,  (2.12) 


where  Qq  is  the  radius  of  the  Earth  and  lO^/cic  =  0.157  m~',  a  is  approximately 
0.120  M-units  per  meter  for  standard  conditions,  d  is  the  duct  height,  and  zq  is  the 


'^T.  E.  Van  Zandt,  J.  L.  Green,  K  S.  Gage,  and  W.  L.  Clark,  "Vertical  Profiles  of 
Refractivity  Turbulence  Structure  Constant-  Comparison  of  Observations  by  the 
Sunset  Radar  with  a  New  Tlieoretical  Model,"  Radio  Sci  13.  819-829  (1978). 

**M.  R  Battaglia,  Modelling  the  Radar  Evaporative  Duct,  Department  of  Defence, 
Defence  Science  and  Teclmology,  RAN  Research  Laboratory,  Australia  (19'i5). 
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"momentum  roughness  length."  Levy  and  Craig  set  a  numerical  value  zq  =  1.5  x 
10"'*  m,  which  is  consistent  with  a  calm  open  sea.  (Panofsky  and  Dutton  have 
tabulated  values  of  zq  for  various  sea  states  and  terrains  [see  Ref.  8,  p.  123].) 
Taking  the  derivative  of  Equation  2.12  and  substituting  the  numerical  constants 
yields 


d<n> 

dz 


'-o.nod 

Z+Zq 


0.037 


10" 


(2.13) 


assuming  a  natural  logarithm  in  Equation  2.12.  Substituting  Equation  2.13  into 
Equation  2.11  gives  the  structure  function.  For  heights  zo«  z  <  d,  zq  and  the 
height-independent  part  of  Equation  2.13  can  be  neglected.  Setting  0  =  =  1 

yields 


Cl{z)  =  a^k^'^{0.n0dfz-^'^10-^^ ,  (2.14) 


which  agrees  with  the  corrected  equation  given  by  Levy  and  Craig.  (Levy  and 
Craig  erroneously  neglect  the  factor  lO”’^  in  their  equation,  but  include  it  in  their 
subsequent  figures.)  Equation  2.14  shows  the  familiar  z"2/3  dependence 
characteristic  of  the  atmospheric  surface  layer. Clearly,  however.  Equation  2.14 
is  valid  only  over  a  limited  regime:  near  the  surface,  zq  becomes  important  to  avoid 
the  singularity  at  z  =  0,  and  for  z  >  d  the  constant  part  of  Equation  2.13  should  be 
included. 

Although  Levy  and  Craig  present  plots  of  the  transmission  loss  versus  height 
and  range  only  up  to  20  m,  it  is  interesting  to  consider  the  behavior  of  the 
structure  constant  at  higher  altitudes  outside  the  surface  area.  Above  the  boundary 
layer,  the  outer  scale  is  approximately  constant  for  z»d  and  Equation  2.9  reduces 
to 


=  (2.  8)lJ'^  (0.037)2 10"'^  (2.15) 


Equating  Equations  2.7  and  2.15  relates  the  variance  of  the  medium  to  the  outer 
scale: 


(n})  =  2.0xl0"*^4- 


(2.16) 


A  typical  value  of  the  variance  is  10~*^,  which  gives  an  outer  scale  of  22  m  and  is 
within  the  range  of  acceptable  values. 
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2.5  LIMITS  OF  SPECTRAL  MODELING 


As  discussed  in  section  2,  the  spectral  modeling  approach  for  generating 
numerical  realizations  of  turbulence  has  a  number  of  attractive  features.  It  has  also 
been  the  traditional  method  used  in  conjunction  with  the  parabolic  equation.^^  It  is 
not,  however,  without  certain  drawbacks  and  limitations.  These  limitations  and 
prospective  alternative  approaches  are  considered  in  the  discussion  that  follows. 

In  the  proposed  spectral  approach  for  generating  realizations,  a  Gaussian  white 
noise  process  is  first  generated  using  a  random  number  generator.  The  resulting 
process  is  then  filtered  in  accordance  with  the  desired  spectrum  to  yield  a  Gaussian 
process.  The  details  of  this  procedure  are  given  in  Appendix  A.  It  is  well  known 
from  probability  theory  that  a  Gaussian  process  is  completely  characterized  by  its 
second  moment  and  associated  spectrum.^*^  Consequently,  when  used  with  the 
parabolic  equation  method,  this  should  produce  reasonable  results  for  the  second- 
order  statistics  of  the  propagating  field. 

Unfortunately,  atmospheric  turbulence  is  typically  an  intermittent  phenomena 
obeying  non-Gaussian  statistics;  areas  of  strong  turbulence  may  be  surrounded  by 
quiescent  regions.  Incorporating  these  features  in  a  model  could  be  crucial  if  the 
objective  of  a  numerical  study  is  to  identify  potentially  catastrophic  events  and  to 
assign  a  probability  to  their  occurrence.  Intermittency  is  primarily  a  factor  at 
higher  altitudes  outside  the  surface  layer,^  but  it  can,  to  a  lesser  extent,  also  occur 
at  lower  altitudes.  Although  filtered  white  noise  models,  as  noted  by  Panofsky  and 
Dutton,  produce  the  correct  average  results,  they  are  incapable  of  simulating  these 
low-probability  events. 

A  possible  alternative  to  spectral  modeling  is  the  so-called  full  simulation 
method  (see  Ref.  9,  pp.  159-161).  In  this  approach,  the  underlying  Navier-Siokes 
equations  that  describe  the  dynamics  of  the  turbulence  are  solved  for  quantities  like 
temperature  and  pressure.  Using  Equation  2.8  or  a  similar  expression,  these 
parameters  are  converted  to  an  index  of  lefr^tion.  The  electromagnetic  propagation 
can  then  be  simulated  using  the  parabolic  equation.  This  approach  was  used  by 
Rouseff  et  al.^^  to  study  the  effects  of  oceanic  microstructure  on  acoustic 
propagation  over  short  ranges.  Although  the  entire  procedure  was  conceptually 
complete,  it  is  likely  to  be  impractical  for  large-scale  atmospheric  simulations. 
The  inner  and  outer  scales  of  turbulence  differ  by  perh^s  five  orders  of  magnitude; 
to  simulate  the  entire  medium  with  the  necessary  resolution  out  to  ranges  on  the 
order  of  tens  of  kilometers  would  be  computationally  prohibitive.  The  increased 
availability  of  super  computing  technology  may  make  this  approach  more 
attractive  in  the  future. 

Wave  propagation  theories  for  random  media  have  traditionally  been  formulated 
in  terms  of  the  spectrum  of  the  medium.  More  recently,  research  has  started  to 


**C.  Macaskill  and  T.  E.  Ewart,  "Computer  Simulation  of  Two-Dimensional  Random 
Wave  Propagation,"  IMA  J.  Appl.  Math.  33,  1-15  (1984). 

^®A.  Papoulis,  Probability,  Random  Variables,  and  Stochastic  Processes,  McGraw- 
Hill,  New  York  (1984). 

Rouseff,  K.  B.  Winters,  and  T.  E.  Ewart,  "Reconstruction  of  Oceanic  Micro¬ 
structure  by  Tomography,"  J.  Geophys.  Res.  96,  8823-8834  (1991). 
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focus  on  calculating  the  probability  distribution  of  the  field.^^  As  these  theories 
become  more  fully  developed,  it  is  possible  that  they  could  be  used  to  assign 
probabilities  of  occurrence  to  abnormal  events. 

A  new  ^rproach  using  fractals  also  has  some  potential  for  modeling  turbulence. 
Kim  and  JaggarcP^  used  fractals  to  simulate  events  exhibiting  both  turbulence-like 
spectra  and  intermittency.  Jaggaid^^  developed  analytical  solutions  for  optical 
propagation  through  a  succession  of  firactal  phase  screens.  Because  of  the  close 
parallels  between  propagation  through  phase  screens  and  what  is  actually 
implemented  in  a  parabolic  equation  routine  (see  section  4.1),  it  seems  likely  that 
fractal  phase  screens  could  be  used  in  TEMPER.  It  appears  unlikely,  however,  that 
fractal  phase  screens  could  be  generated  using  efficient  numerical  techniques  like 
the  FFT. 

To  summarize,  spectral  models  are  relatively  easy  to  implement  and  have 
traditionally  been  used  in  conjunction  with  the  parabolic  equation  method.  Useful 
results  have  long  been  obtained  in  ocean  acoustics  rqrplications.  Spectral  modeling 
can  be  a  fruitful  approach  if  the  objective  of  a  numerical  study  is  to  calculate  the 
average  fluctuations  in  the  field  caused  by  randomness.  The  method  cannot, 
however,  be  used  to  simulate  low-probability  events. 


E.  Ewart,  "A  Model  of  the  Intensity  Probability  Distribution  for  Wave  Propa¬ 
gation  in  Random  Media,”  J.  Acoust.  Soc.  Am.  86,  1490-1498  (1989). 

^^E.  Jakeman,  "Scattering  by  Gamma-Distributed  Phase  Screens,"  Waves  in  Random 
Media  2,  153-167  (1991). 

Kim  and  D,  L.  laggard,  "Band-Limited  Fractal  Model  of  Atmospheric  Re- 
fractivity  Fluctuation,"  J.  Opt.  Soc.  Am.  A:  Opt.  Image  Sci.  5,  475-481  (1988). 

2*D.  L.  laggard,  "On  Fracul  Electrodynamics,"  in  Recent  Advances  in  Electro¬ 
magnetic  Theory,  H.  N.  Kritikos  and  D.  L.  laggard  (ed.),  Springer-Verlag,  New 
York  (1991). 
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3.  THEORETICAL  SOLUTIONS  FOR  PROPAGATION 
THROUGH  RANDOM  MEDIA 

Tractable  analytical  solutions  for  predicting  the  statistical  properties  of  an 
electromagnetic  wave  propagating  through  a  random  medium  are  available  only  for 
certain  idealized  source  configurations  such  as  a  plane  wave  or  a  point  source. 
These  solutions  usually  neglect  both  deterministic  index  of  refraction  profiles  and 
the  interaction  of  the  field  with  any  boundaries.  To  solve  these  more  difficult 
problems,  numerical  techniques  such  as  the  paraljolic  equation  method  are 
necessary.  For  two  important  reasons,  however,  it  is  still  cruc'al  to  consider  the 
analytical  solutions,  despite  their  limited  applicability.  First,  the  solutions  for  the 
canonical  problems  provide  a  benchmark  for  testing  numerical  algorithms. 
Second,  the  analytical  solutions  provide  insight  as  to  how  three-dimensional 
turbulent  fluctuations  should  be  modeled  for  inclusion  in  two-dimensional 
numerical  simulations. 

Fundamentally,  TEMPER  solves  a  two-dimensional  wave  propagation 
problem.  This  two-dimensional  problem  is  interpreted  as  representing  a  cross- 
sectional  slice  of  the  true  three-dimensional  reality,  which  is  reasonable  for 
propagation  through  deterministic  media  that  are  varying  in  a  scale  very  large 
compared  to  the  wavelength.  Random  turbulent  fluctuations,  however,  can  vary  on 
scales  as  small  as  a  millimeter.  Consequently,  the  relationship  between  the  two- 
and  three-dimensional  problems  is  no  longer  clear.  In  section  3.4,  this  relationship 
is  examined  in  detail  for  plane  wave  propagation.  If  the  statistical  quantity  of 
interest  is  the  second  moment  of  the  field,  then  the  two-dimensional  result 
accurately  predicts  the  cross-sectional  behavior  of  the  true  three-dimensional 
problem.  Two-dimensional  simulations  cannot,  however,  accurately  predict  the 
true  log-amplitude  fluctuations  except  in  the  very  far  field. 

3.1  STRONG  AND  WEAK  SCATTERING  THEORIES 

First,  neglect  atmospheric  ducting  and  sea  surface  interaction  and  consider  the 
idealized  problem  of  a  wave  propagating  through  a  medium  characterized  by  small 
random  fluctuations  in  its  index  of  refraction.  After  propagating  a  short  distance, 
the  wavefront  will  start  to  acquire  small  deviations  from  its  expected  value  owing 
to  the  randomness.  The  total  field  can  be  interpreted  as  the  original  coherent  wave 
front  plus  a  small  random  incoherent  component.  At  greater  distances,  the 
cumulative  effect  of  the  random  scattering  will  produce  fluctuations  on  the  same 
order  as  the  original  signal.  Finally,  at  still  greater  ranges,  the  random  component 
completely  dominates  and  the  field  is  essentially  incoherent. 

The  theoretical  approach  used  to  calculate  the  statistical  properties  of  the 
propagating  wave  depends  on  the  range  of  interest.  At  sufficiently  short  ranges, 
where  the  contribution  due  to  randomness  is  small,  weak  scattering  theory  is 
applicable.  Weak  scattering  theory  uses  either  the  Born  or  the  Rytov 
approximation  to  model  the  propagating  field  within  the  medium.  For  forward 
propagation  through  smoothly  fluctuating  media,  the  Rytov  approximation  is 
generally  regarded  as  superior.  If  the  scattering  is  sufficiently  weak,  there  is  no 
practical  difference  between  the  two  approximations.  At  longer  ranges,  where  the 
incoherent  component  of  the  field  is  large,  strong  scattering  theory  is  required.  In 
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this  approach,  partial  differential  equations  are  derived  to  describe  the  first,  second, 
and  higher  moments  of  the  total  field.  Beyond  the  second  moment,  numerical 
techniques  are  generally  necessary  to  solve  the  equations  approximately;  this  is  the 
subject  of  active  research. 

Microwave  propagation  through  a  random  medium  differs  from  the  analogous 
optical  and  ocean  acoustic  scattering  problems  in  that  most  practical  problems  can 
be  modeled  by  weak  scattering  theory..  Figure  3.1  shows  the  rough  regions  of 
validity  of  the  alternate  theoretical  approaches  for  plane  wave  propagation  through 
a  strongly  turbulent  atmosphere.  At  optical  frequencies,  strong  scattering  theory  is 
necessary  beyond  about  1  km..  At  a  frequency  of  10  GHz,  however,  weak 
scattering  is  valid  to  a  height  of  roughly  1000  km.  Consequently,  weak  scattering 
theory  will  be  emphasized  in  this  analysis.  It  can  be  shown  that  strong  scattering 
theory  produces  expressions  that  are  consistent  with  weak  scattering  predictions  in 
the  appropriate  limiting  cases  and  for  the  appropriate  statistical  quantities. 

Weak  scattering  theory  gives  integral  expressions  for  various  statistics  of  the 
field  in  terms  of  the  fluctuations  of  the  medium.  The  approach  used  to  solve  the 
integrals  depends  on  the  frequency  and  range  of  interest,  and  also  on  the  scale  of  the 
irregularities.  As  noted  in  section  2,  turbulence  is  characterized  by  an  inner  scale 
^0  and  an  outer  scale  Lq-  Figure  3.1  shows  the  subregions  of  weak  scattering  for 
various  values  of  Iq  and  Lq.  Consider  a  range  L  and  an  illuminating  wavelength 
A.  If  L  «  /  A,  the  Rytov  integral  is  approximated  using  geometric  optics;  for 

microwave  frequencies,  this  region  is  of  little  practical  interest.  For  ranges 
^0  /  A  «  L  «  Zfl  /  A,  both  the  inner  and  outer  scales  can  affect  the  propagation 
and  the  field  is  said  to  be  in  the  diffraction  (Fresnel)  regime.  (A  note  of  caution  is 
in  order.  The  phrase  "diffraction  regime"  is  used  here  differently  than  in  radar  work, 
where  it  commonly  refers  to  over-the-horizon  propagation.')  Finally,  for 
L»f^/A,  the  Rytov  integral  is  well  approximated  using  the  far  field  or 
Fraunhofer  approximation. 

Let  the  field  at  point  r  after  propagating  from  the  source  through  randomness  be 
given  by  U{v)  =  A(r)  exp  [(/Fr)].  Designate  the  field  that  would  exist  at  the  same 
point,  coming  from  the  same  source  but  without  randomness,  as  Voir)  = 
Aq  exp  [j7’o(r)]'  As  noted  earlier,  weak  scattering  theory  gives  expressions  for 
various  statistical  quantities.  Among  these  are  the  autocorrelation  of  the 
normalized  log-amplitude  of  the  field  and  the  autocorrelation  Bp  of  the  normalized 
phase: 


S;t(n.r2)  =  (A(n)A(r2)) 

(3.1a) 

Bp(r,,r2)  =  (P](rj)F,(r2)), 

(3.1b) 

where  the  angle  brackets  indicate  an  ensemble  average  and 
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Figure  3.1  Scattering  regimes  for  a  piane  wave  propagating  through  strong 
turbuience.  The  negativeiy  sioping  shaded  iine  indicates  the  break  between  strong 
and  weak  scattering  theories.  Weak  scattering,  based  on  the  Rytov 
approximation,  is  further  divided  into  the  geometric  optics,  the  diffraction,  and  the 
Fraunhofer  subregions.  The  broken  iines  between  the  weak  scattering  subregions 
depend  on  the  inner  scale  /q  4  of  the  turbulence.  The  line 

between  strong  and  weak  scattering  is  set  where  the  variance  of  the  field  log- 
amplitude  fluctuations  equals  0.2  (see  Eq.  3.9)  and  the  structure  constant  for 
turbulence  is  set  to  1 0”^. 
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;i;(r)  =  ln[A(r)/Ao(r)] 


(3.2a) 


Pi(r)  =  P(r)-Po(r). 


(3.2b) 


The  log-amplitude  is  frequently  a  quantity  of  interest  for  radar  work.  Since  the 
phase  is  related  to  the  travel  time  of  the  wave,  the  phase  fluctuations  are  often  the 
quantities  of  interest  in  ocean  acoustics.  The  associated  structure  functions  are 
given  by 


Dx(ri,T2)  =  (jx0ri)-Xir2f)  (3.3a) 

Dp(ruT2)  =  {\f\(ri)-Pi(T2f).  (3.3b) 


Both  the  log-amplitude  and  the  phase  are  quantities  that  are  natural  to  calculate 
within  the  Rytov  approximation.  Within  the  Bom  approximation,  ic  is  natural  to 
consider  the  scattered  field 


l/,(r)=£/(r)-[/o(r)  (3.4) 

with  the  associated  autocorrelation 

BM>r2)  =  {UsWUr2)).  (3.5) 


Although  the  log-amplitude  and  phase  are  both  real  functions  of  position,  the 
scattered  field  is  complex.  Note  also  that  the  mean  incident  field  is  subtracted  from 
the  total  in  Equation  3.4  to  yield  the  scattered  field.  Hence  for  weak  scattering  the 
autocorrelation  of  the  scattered  field  can  be  interpreted  as  the  autocovariance  of  the 
total. 

Strong  scattering  theory  produces  expressions  for  the  statistical  moments  of  the 
total  complex  field  f/(r)  =  A(r)  exp(<P(r)].  Much  of  the  current  research 
concentrates  on  solving  for  the  fourth  moment,  which  can  be  related  'o  intensity 
scintillation.  For  comparison  with  weak  scattering  theory,  however,  the  second 
moment  of  the  field  will  be  emphasized  at  range  L,  defined  by 


r(x  =  L\yi,2i’,y2,zz)  =  {U{x  =  L,yi,zi)U*{x  =  L,y2,Z2)).  (3.5) 


Particularly  in  optics,  the  second  moment  is  called  the  muUtal  coherence  function. 
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3.2  PLANE  WAVE  PROPAGATION 


Consider  the  correlation  of  the  field  at  points  rj  and  r2  at  common  range  x\  = 
JC2  =  L.  If  the  medium  is  statistically  homogeneous,  the  field  statistics  will  be  a 
function  of  the  difference  coordinates  yd= yi-  yi  and  zj  =  zj  -  22-  The  integral 
expressions  within  the  Rytov  approximation  are^® 


5;^ (f-. )'<i . dk^ S„ (Q,ky,k^ ) cxip[-i(kyyj  +  k^zj )]f^ ( k) , 

(3.7a) 


Bp{L,ya,zj)  =  7[k^Lf  dkyf  dk2S„(0,ky,ki)c\p[-i{kyyd  + k,Zd)]fp{K), 

(3.7b) 


0  'y 

where  k  -InlXis  the  wave  number,  k  -ky+k^,  and  S„  is  the  three- 
dimensional  spectrum  of  the  index  of  refraction.  The  functions and  fp  are  called 
the  spectral  filter  functions  for  the  log-amplitude  and  phase,  respectively,  and  are 
given  by 


/,(.)=  1-^^,  (3.8a) 

=  (3.8b) 


If  the  medium  is  also  isotropic,  single  integr'd  expressions  are  available  (see  Ref., 
4,  p.  359). 

For  a  turbulent  medium  within  the  Fresnel  regime,  tl/  X«L«L^/  X,  the 
dominant  contribution  to  the  integration  in  Equation  3.7a  comes  from  the  inertial 
subrange  of  the  Kolmogorov  spectrum  and  S„  =  0.033)f’''^  can  be  used  for  all  k. 
Evaluating  at  =  Zd=  0  gives  the  variance  of  the  log-amplitude  (see  Ref.  4,  p. 
370): 


ol  =  B^{L,yj  =  Zd  =  0)  =  O.^OlClk'^'^L^  .  (3.9) 


Equation  3.9  was  used  to  define  the  transition  between  weak  and  strong  scattering 
in  Figure  3.1.  Weak  scattering  theory  can  be  used  for  <  0. 2  -  0. 5 .  Within 
the  Fresnel  regime.  Equation  3.7a  can  also  be  approximated  for  arbitrary  correlation 
distance.  The  resulting  expression  is  exU’emeiy  cumbersome  (see  Ref.  4,  p.  370). 

The  general  expressions  for  the  autocorrelations  simplify  considerably  within 
the  Fraunhofer  regime.  Except  for  very  small  k,  the  spectral  filter  functions  of 


M.  Rytov,  Y  A.  Kravtsov,  and  V  1.  Tatarskii,  Principles  of  Statistical 
Radiophysics,  Vol.  4,  Springer-Verlag,  New  York  (1989). 
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Equation  3.7  are  both  approximately  equal  to  I  for  the  significant  part  of  the 
integration  (see  Ref.  4,  p.  364,  and  Ref.  26,  p.  57): 


f^=fp=l  (Fraunhofer regime).  (3.10) 

Equation  3.7  can  thus  be  approximated  by 

=  Bp  =k^L^  B„(Xd,yd,2d)dXd,  (3.11) 

where  it  has  been  assumed  that  the  autocorrelation  of  the  medium  exists. 

Approximating  the  spectral  filters  by  a  constant  simplifies  the  resulting  integral 
expressions  for  the  log-amplitude  and  phase  fluctuations.  However,  the 
approximate  expressions  are  only  marginally  useful  for  propagation  through 
turbulence.  For  example,  consider  a  3.3-GHz  plane  wave  propagating  through  a 
turbulent  medium  with  an  outer  scale  of  100  m.  The  Fraunhofer  criterion  is  met 
for  ranges  much  greater  than  110  km.  For  an  outer  scale  of  10  m,  the  far  field 
applies  at  ranges  much  greater  than  1.1  km. 

For  the  scattered  field,  the  autocorrelation  is  given  by 

BsiL,yd,2d)  =  nk^Lf‘  dkyf  dktS„iO,ky,ki)Qxi[>[-i(kyyd  + 

(3.12a) 

where  /„  the  spectral  filter  for  the  scattered  field,  can  be  writteii  in  terms  of  the 
specU'al  filters  for  the  log-amplitude  and  the  phase: 


fs=fx+fp  =  2.  (3.12b) 

Note  that  the  specU'al  filler  for  the  scattered  field  is  constant  fcr  arbitrary  range. 
This  is  in  contrast  to  the  log-amplitude  and  the  phase,  where  the  the  associated 
spectral  filters  are  approximately  constant  only  at  very  distant  ranges.  Assuming 
the  autocorrelation  of  the  medium  exists.  Equation  3.12a  can  be  rewritten  as 

Bs{L,yd,Zd)  =  2k^L^B„  (x^ .  )£iCd .  (3.13) 

Equation  3.13  will  be  used  in  section  4  to  calculate  the  expected  solutions  to 
which  numerical  simulations  will  be  compared. 

For  strong  scattering  theory,  the  mutual  coherence  function  is  given  by  (see 
Ref.  4,  p.  414,  and  Ref.  6,  p.  31) 


in=Qexp-i£>(>^,zj)j 


(3.14) 


where  =  1  for  a  plane  wave.  Diy^.zd)  is  the  wave  structure  function  and  is 
related  to  the  log-amplitude  and  phase  structure  functions  used  in  Rytov  theory  (Eq. 
3.2)  by 


.  z<i)  =  %  (yrf  .Zrf ) + Dp  {yd ,  z<i ).,  (3.15) 


By  appropriate  expansions,  it  can  be  shown  than  Equation  3.14  is  consistent  with 
Equation  3.13  when  IDI «  1. 

Equation  3.14  is  a  typical  strong  scattering  result  When  the  random  fluctua¬ 
tions  in  the  index  of  refraction  are  assumed  to  occur  about  some  deterministic  mean 
profile,  the  wave  structure  function  can  be  written  in  terms  of  path  integrals 
through  the  medium.^^'^*  Evaluation  of  these  path  integrals  is  the  subject  of 
current  research  in  the  scattering  community. 

3.3  POINT  SOURCE  PROPAGATION 

Point  source  illumination  is  another  important  geomeuy  for  which  closed  form 
expressions  are  available.  Because  the  numerical  simulations  in  section  4 
emphasize  plane  waves,  only  the  primary  point  source  results  are  summarized  here. 
(For  greater  detail,  see  Refs.  4  and  7.)  In  weak  scattering,  the  derivation  begins  by 
making  the  Rytov  approximation  to  the  scattering  integral.  A  paraxial  approxima¬ 
tion  is  then  made  for  both  the  incident  field  and  the  free  space  Green’s  function. 
The  variance  of  the  normalized  log-amplitude  fluctuations  for  a  turbulent  medium 
is  given  by  (Ref.  4,  p.  379) 


(3.16) 


which  can  be  compared  to  the  analogous  plane  wave  result  in  Equation  3.9. 

For  suong  scattering,  the  mutual  coherence  is  again  given  by  Equation  3.14 
with  Cs  =  L~'^. 

3.4  RELATIONSHIP  BETWEEN  TWO-  AND  THREE- 
DIMENSIONAL  PROBLEMS 

The  parabolic  equation/split-step  algorithm  as  implemented  in  TEMPER  solves 
a  two-dimensional  propagation  problem.  The  results  are  interpreted  as  representing 
propagation  in  a  vertical  cross-section  of  the  true  three-dimensional  geometry. 


M.  Flalte,  R.  Dashen,  W.  H.  Munk,  K.  M.  Watson,  and  F.  Zachanasen,  Sound 
Transmission  Through  a  Fluctuating  Ocean.  Cambridge  University  Press,  New  York 
(1979). 

M.  Flatte,  "Wave  Propagation  through  Random  Media'  Contributions  from 
Ocean  Acoustics,"  Proc  IEEE  71,  1267-1294  (1983). 
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This  is  reasonable  for  propagation  through  a  deterministic  medium,  since  it  is 
typically  fluctuating  on  a  scale  that  is  very  large  compared  to  the  wavelength  A  and 
out-of-plane  scatter  is  small.  It  is  not  obvious,  however,  that  this  approximation 
remains  valid  when  random  turbulent  fluctuations  are  included.  As  noted  in  section 
2.1,  the  inner  scale  of  turbulence  can  be  as  little  as  a  millimeter.  Consequently,  a 
turbulent  medium  will  contain  fluctuations  that  are  small  compared  to  A  at 
microwave  frequencies.  In  this  section,  theoretical  solutions  for  the  two- 
dimensional  problem  that  is  actually  simulated  by  TEMPER  are  considered.  The 
conditions  under  which  the  two-dimensional  solutions  might  be  legitimately 
interpreted  as  representing  a  cross-section  of  the  three-dimensional  problem  are 
given. 

Consider  the  statistics  of  the  propagating  field  evaluated  along  a  vertical  line. 
The  weak  scattering  plane  wave  results  given  in  Equations  3.7  and  3.12  can  be 
written  in  a  compact  notation  as 


BQ(L,yj,2^)  =  7tk^L  dky 


dkX{Q,ky,k^)  t\'p[-i{kyyd  +  k^Zd)]fQ{K) , 

(3.17) 


where  g  can  be  ;i;,  P,  or  s  to  represent  the  log-amplitude,  the  phase,  or  the 
scattered  field,  respectively.  Following  the  derivation  of  Ishimaru,'*  it  can  be 
shown  that  the  same  statistical  quantities  in  the  two-dimensional  problem  are 
given  by 


B^'>{x  =  L,2d)  =  nk 


h 


dk^Sj^\k^  =  O.ki)  t\p{+ikiZ^)fQ{ki).-  (3.18) 


Comparing  the  two  expressions,  we  see  that  Equation  3.17  reduces  to  Equation 
3.18  for  any  specual  filter  function  /g  if  the  three-dimensional  spectrum  is  of  the 
form 


S,{k,,ky,k,)=Sj^\kMS{ky).  (3.19) 


Equation  3.19  states  the  trivial  result  that  the  three-dimensional  problem  reduces  to 
two  dimensions  if  there  is  no  variation  in  the  transverse  coordinate:  this  clearly  is 
not  the  case  for  atmospheric  turbulence.  As  an  alternative,  we  must  design  the 
two-dimensional  spectrum  to  mimic  the  three-dimensional  fluctuations.  The  way 
in  which  this  is  done  will  depend  on  the  spectral  filter  function /g.  For  the 
scattered  field,  g  =  ^  and  the  spectral  filter  is  constant.  Under  these  circumstances, 
Equations  3.17  and  3.18  are  equivalent  a.tya  =  0  if 

sl^\k„k,)=  r  dkyS„{k„ky,k,).  (3.20) 

J  —00 


Taking  the  inverse  uansform  of  Equation  3.20  relates  the  two  autocorrelations 
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=  Bn(Xd>yd  = 


(3.21) 


Hence  the  two-dimensional  random  process  is  simply  a  cross-section  of  the  three- 
dimensional  process  evaluated  along  a  vertical  plane. 

Equation  3.21  is  a  significant  result;  it  implies  that  meaningful  results  can  be 
obtained  from  two-dimensional  simulations.  Further,  it  gives  a  simple 
relationship  between  the  three-dimensional  statistics  of  the  medium  that  are 
measured  and  the  two-dimensional  statistics  that  will  be  used  in  the  simulations. 
This  cross-sectional  notion  has  been  routinely  used  in  numerical  studies  of  strong 
scattering  phenomena.'^  This  analysis  suggests  that  this  assumption  is  also  true 
for  weak  scattering  provided  the  statistical  quantity  of  interest  is  the  fluctuations  of 
the  scattered  field. 

Mathematically,  this  dimensional  simplification  was  possible  because  the 
specU'al  filter  for  the  scattered  field  is  a  constant.  Except  in  the  far  field,  however, 
the  spectral  filter  functions  for  the  log-amplitude  and  the  phase  fluctuations  are  not 
constant  and  the  simple  partitioning  used  above  is  not  valid.  At  ranges  in  the 
Fresnel  regime,  the  oscillatory  part  of  the  spectral  filters  in  Equation  3.8 
contributes  to  the  autocorrelation  integrals.  At  distant  ranges  in  the  Fraunhofer 
regime,  the  oscillatory  contribution  can  be  neglected,  the  spectral  filters 
approximated  by  constants  (Eq.  3.10),  and  the  cross-sectional  assumption  again 
becomes  reasonable. 

To  summarize,  if  the  statistical  quantity  of  interest  is  the  autocorrelation  of  the 
scattered  field,  then  it  is  valid  to  simulate  a  cross-section  of  the  three-dimensional 
medium.  The  result  holds  for  arbitrary  range.  In  calculations  of  the  log-amplitude 
fluctuations  or  the  phase  fluctuations,  however,  the  cross-sectional  assumption  is 
approximately  correct  only  in  the  far  field. 

To  quantify  the  results  of  the  preceding  section,  consider  a  medium  characterized 
by  an  isotropic  Gaussian  autocorrelation: 


Bn{x^,yd’^d)=(^l  exp[-(;cj  +  +  Zd )  /  f  o  ]  •  (3.22) 


While  the  Gaussian  spectrum  does  not  represent  any  uue  medium,  it  is 
mathematically  more  convenient  to  manipulate  than  the  von  Kiuman  spectrum  and 
is  adequate  for  this  example.  An  expression  for  the  variance  of  the  log-amplitude 
is  given  by  Ref.  26  (p.  77).  Converung  to  the  notation  used  in  this  paper. 


ol=^^al  k^L fo[l  -  £>■'  tan-'  (D)] ,  (3.23) 

where  D  =  kLI{k(Ql2)  .  The  first  term  inside  the  square  brackets  in  Equation 
3.23  is  from  the  constant  part  of  the  spccual  lilter  and  is  the  far  field  (large  D) 
result.  The  second  term,  D"'  ian~\D),  is  due  to  the  oscillatory  part  of  the  spectral 
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filter;  it  can  be  interpreted  as  the  Fresnel  (diffraction)  correction.  For  ranges 
L »  /§  /  ^ .  that  is,  D»  I,  this  term  falls  like  £>"*. 

A  similar  derivation  is  possible  for  a  two-dimensional  geometry.  Evaluating 
Equation  3.22  at  =  0  yields  the  two-dimensional  autocorrelation.  Taking  the 
Fourier  transform  yields  the  spectrum.  The  spectrum  along  with  the  log-amplitude 
spectral  filter  (Eq.  3.8a)  are  substituted  into  Equation  3.18.  With  =  0  and  a 
suitable  change  of  variables,  the  resulting  integral  can  be  evaluated,^^  yielding 


2D~\l  +  D^)‘/'*  sin 


^tan"* 


(!»]}. 


(3.24) 


As  expected,  the  first  term,  that  is,  the  fiu*  field  Fraunhofer  result,  is  equivalent  in 
the  two  expressions.  Note,  however,  that  the  second  term  falls  like  for 
D  »  1  This  is  a  significantly  slower  falloff  than  is  observed  for  the  three- 
dimensional  case;  the  two-dimensional  case  reaches  the  far  field  more  slowly. 

As  a  numerical  example,  consider  the  following  typical  values: 


Frequency 
Range 
Correlation  length 
Variance  of  medium 
D 


3.3  GHz  (t  =  69.12  m->) 

4.9  nautical  miles  (L  =  9075  m) 
10  feet  (/o  =  3.04  m) 
cr2  =  10-‘2 

()tL)/(/t  4/2)2  =  56.8 


Substituting  into  Equations  3.23  and  3.24  yields 


Three-dimensions: 

^  Vw  <T^  jfe  fo  (1  -  0- 027)  =  1 1. 36  X 1 0*^ 


Two-dimensions: 

oj  =^^niolk^Lto{\-0.m)  =  9.5lxlO-\ 

While  the  correction  term  in  the  three-dimensional  case  makes  only  a  2.7% 
contribution,  the  corresponding  term  in  two-dimensions  contributes  18.6%.  This 
example  suggests  that  the  use  of  two-dimensional  simulations  to  calculate  the 
three-dimensional  log-amplitude  fluctuations  must  be  viewed  skeptically. 


S.  Gradshteyn  and  1.  M.  Ryzhick,  Tables  of  Integrals,  Senes  and  Products, 
Academic  Press,  Orlando,  Fla.,  p.  490  (1980). 
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4.  MODELING  TURBULENT  FLUCTUATIONS 
IN  TEMPER 


Having  developed  the  necessary  background  theory,  we  are  now  ready  to 
consider  proposed  methods  for  including  turbulent  refractive  index  fluctuation  in 
TEMPER. 

As  shown  by  Ko  et  al.,^®  derivation  of  the  electromagnetic  parabolic  equation 
usually  begins  by  assuming  a  spherical  coordinate  system  with  the  origin  at  the 
center  of  an  Earth  of  radius  ag  and  a  vertical  electric  (or  magnetic)  dipole  source 
located  at  r  =  a,  +/i  and  0=  0.  In  deriving  the  time-independent  wave  equation,  all 
derivatives  with  respect  to  the  0  coordinate  are  ignored.  This  is  justified  through 
the  rotational  symmetry  of  the  source  and  by  assuming  the  medium  is  rotationally 
symmetric  about  the  source  point.  The  problem  is  reduced  to  two  dimensions  and 
through  a  conformal  transformation  it  is  written  in  pseudo-rectangular  coordinates.^ 
Some  form  of  the  parabolic  approximation  is  then  made  to  admit  a  marching 
solution  through  the  Fourier  split-step  method.  Using  this  two-dimensional 
assumption  as  a  starting  point,  a  method  for  including  random  fluctuations  is 
developed. 

4.1  INCLUDING  RANDOMNESS  IN  TEMPER 

TEMPER  solves  the  two-dimensional  parabolic  equation  with  height  z  and 
range  x.  The  field  u  at  range  x  is  propagated  to  x  +  &  by  the  algorithm' 


u(x  +  5x,  z)  =  ^  exp[i/:0(z)] 


r°° 


dp  exp(i/jz)  exp 

✓  — oo 


-ip^Sx 

~ir~ 


dz'u{x,z')  exp(-/pz'), 

J  —00 


(4.1) 


where  k  =  2n/X  is  the  free  space  wave  number  and  exp[/I:0(z)]  is  the  transmittance 
function  with 


,  rX-*-Sx  .  , 

0(z)  =  -jJ  dx’^n^{x',z)-lj  + z6x  /  Qg ,  (4-2) 

where  n  is  the  index  of  refraction. 

The  procedure  can  be  summarized  as  follows:  lake  the  Fourier  transform  with 
respect  to  height  of  the  given  initial  field  distribution,  multiply  by  a  propagauon 
filter,  inverse  nansform,  and  finally  multiply  by  a  transmittance  function  that 
represents  the  effects  of  the  medium  between  x  and  .r  +  5x.  The  split-step 
algorithm  separates  the  propagation  effects  from  the  effects  of  the  medium;  the 
field  is  propagated  a  distance  5x  as  if  through  free  space  and  then  a  posteriori 
multiplied  by  the  u-ansmiitance  that  modulates  the  phase  of  the  field  to  account  for 


W.  Ko,  J.  W.  Sari,  and  J.  P.  Skura,  "Anomalous  Microwave  Propagation 
through  Atmospheric  Ducts,"  Johns  Hopkins  APL  Tech  Dig.  4,  12-26  (1983) 
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the  integrated  effects  of  the  medium.  In  analogy  with  optics,  the  transmittance 
function  is  usually  called  a  phase  screen. 

The  index  of  refraction  includes  a  deterministic,  height-dependent  mean  profile 
<n(z)>  and  a  small  randomly  fluctuation  pannf{x,  z); 


n{x,  z)  =  1+  <  n(z)  >  +ny  (x,  z) , 


(4.3) 


The  profile  may  be  weakly  range  dependent,  but  for  simplicity  is  assumed  here  to 
be  independent  of  x  over  the  step  size.  Assuming  l<«>  +  «  1  and  neglecting 

second-order  terms  yields 


0(z)  =  [<rt>-fz/fl^]&  +  0,(z),  (4.4a) 


where  <t>,  represents  the  integrated  random  refractive  index  fluctuations. 


dx'nf(x\z). 


(4.4b) 


The  objective  is  to  generate  random  realizations  of  0,  consistent  with  the 
assumed  statistics  of  n/.  For  this  two-dimensional  problem,  the  random  part  of  the 
refractive  index  is  modeled  as  a  zero  mean,  statistically  homogeneous  process  with 
an  autocorrelation  function  It  follows  from  Equation  4.4b  that  0,  is 

zero  mean,  with  autocorrelation  given  by 

/’X+&C  rx-^dx 

(0,(zi)0,(z2))  =  J  dxi  J  (ix2fi^^^(xi-X2,r, -Zz).  (4-5) 

Changing  the  integration  variables  to  Xs  =  (xj  +  xi)!!  and  x^  =  (xj  -  xi),  and 
assuming  that  the  step  size  8x  is  large  compared  to  the  correlation  length  of  the 
medium,  yields  approximately 


(0/(zi)0/(z2))  =  5-<r  dXdB^^\xd,Zd)sBt{zd). 

J  — oo 


(4.6> 


The  quantity  B,  is  called  the  uansverse  autocorrelation  and  plays  an  important  role 
in  strong  scattering  theory.*  With  an  eye  toward  eventual  implementation  in 
Appendix  A  using  the  fast  Fourier  U'ansform,  we  define  the  associated  transverse 
spectrum  by 


S,{k,)  = 


dZdS;(Zd)  GXpi-ik^Zd). 


(4.7) 
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For  consistency  with  the  three-dimensional  spectra  used  in  section  2,  the  two- 
dimensional  spectrum  of  the  medium  is  defined  by 


2 

J  — OA  J  — oo 


(4.8) 


Combining  Equations  4.6-4.8  relates  the  transverse  spectrum  to  the  two- 
dimensional  spectrum  of  the  medium 


S,{k,)  =  {2nfSxS^^\k,=0,k,).  (4.9) 


Realizations  of  <>/  consistent  with  this  transverse  spectrum  can  be  generated  using 
the  algorithm  in  Appendix  A.  These  fluctuations  can  then  be  included  in  the 
marching  algorithm  of  Equation  4.1  to  simulate  propagation  through  a  single 
realization  of  a  random  medium.  The  statistical  properties  of  the  propagating  field 
for  any  given  realization  can  be  estimated.  As  part  of  a  Monte  Carlo  procedure,  the 
simulation  can  be  repeated  for  many  independent  realizations.  By  averaging  the 
estimates,  an  overall  estimate  of  the  true  ensemble  characteristics  is  obtained. 

One  of  the  central  assumptions  of  this  derivation  is  that  the  step  size  Sx  is 
large  compared  to  the  outer  scale  of  the  medium  Lq.  With  this  assumption,  the 
phase  screens  at  adjacent  range  steps  are  essentially  unco^elated  and  can  be 
generated  independently.  If  Sx  <  Lq,  then  the  individual  screens  are  correlated  and 
the  medium  must  be  generated  using  a  two-dimensional  FFT  routine.  Spivak^* 
has  shown  that  if  the  objective  of  a  numerical  study  is  to  predict  the  correct 
ensemble  statistics,  then  uncorrelated  phase  screens  are  adequate.  Correlated  screens 
are  necessary  if  the  objective  is  to  mimic  the  fine  detail  that  would  truly  exist  in 
individual  realizations.  Correlated  screens  were  also  used  by  Rouseff  and  Porter^^ 
to  test  stochastic  inverse  scattering  algorithms.  Since  the  generation  of  correlated 
phase  screens  is  much  more  computationally  expensive,  and  since  using  filtered 
white  noise  to  model  turbulence  is  at  best  valid  only  in  an  ensemble  sense  anyway 
(section  2.5),  uncorrelated  screens  will  be  used  in  this  study. 

The  basic  split-step  algorithm  of  Equation  4. 1  can  be  extended  to  three  dimen¬ 
sions.  The  one-dimensional  forward  and  inverse  transforms  are  replaced  by  two- 
dimensional  bansforms  along  the  transverse  y-z  plane.  The  dimensionality  of  the 
phase  screens  is  also  increased.  This  approach  was  used  by  Martin  and  Flatt^^^  to 
model  optical  propagation  through  turbulence.  Unfortunately,  adding  the  extra 
dimension  makes  the  problem  very  computationally  intensive:  Martin  and  Flatte 
required  a  Cray  super  computer  to  implement  their  algorithm.  As  computer 


^*M.  Spivack,  "Accuracy  of  the  Moments  from  Simulation  of  Waves  in  Random 
Media,”  J.  Opt  Soc  Am  A  1 790-793  (1990). 

Rouseff  and  R.  P.  Porter,  "Diffraction  Tomography  and  the  Stochastic  Inverse 
Scattering  Problem,"  J  Acoust  Soc.  Am.  89,  1599-1605  (1991) 

M.  Martin  and  S.  M  Flatte,  "Intensity  Images  and  Staustics  from  Numerical 
Simulation  of  Wave  Propagation  in  3-D  Random  Media,"  Appl  Opt.  21,  2111- 
2127  (1988) 
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technology  improves,  full  three-dimensional  simulations  may  become  practical  for 
desk-top  calculations.  In  the  meantime,  we  are  left  with  two-dimensional 
numerical  models  that  attempt  to  predict  three-dimensional  phenomena. 

The  relationship  between  the  two-  and  three-dimensional  problems  was 
discussed  in  section  3.4.  If  the  quantity  of  interest  is  the  scattered  field,  then  at 
least  for  plane  wave  propagation  it  suffices  to  simulate  a  two-dimensional  cross- 
section  of  the  true  three-dimensional  medium.  If  the  statistical  quantity  of  interest 
is  the  log-amplitude  or  phase  fluctuations  of  the  field,  then  other  strategies  for 
simulating  the  medium  must  be  developed.  For  simplicity,  this  preliminary 
numnical  study  will  concentrate  on  calculating  the  second  moment  The  two-  and 
three-dimensional  spectra  for  the  medium  are  related  by  Equation  3.20: 

S^^^(kM  =  f_^dkyS,(k„ky,k,). 


Combining  Equations  3.20  and  4.9  relates  the  one-dimensional  transverse  spectrum 
used  in  parabolic  equation  simulations  to  the  three-dimensional  spectrum  of  the 
medium: 


S,{k,)  =  {2nfSx£^dkySnik,  =  0,ky,k,).  (4.10) 

4.2  NUMERICAL  EXAMPLE:  PROPAGATION  THROUGH 
A  GAUSSIAN  SPECTRUM 

As  a  numerical  example,  we  consider  plane  wave  propagation  through  a 
medium  described  by  a  Gaussian  spectrum.  Although  the  Gaussian  spectrum  does 
not  describe  an  actual  medium,  its  properties  lead  to  both  a  simple  expression  for 
the  transverse  spectrum  and  relatively  rapid  convergence  in  the  Monte  Carlo 
simulations. 

The  autocorrelation  for  an  isotropic  Gaussian  correlated  medium  is  given  by 
Equation  3.22.  The  ccmesponding  three-dimensional  spectrum  is  easily  calculated 
and  substimted  into  Equation  4.10  to  yield  the  transverse  spectrum. 

In  the  simulations,  a  3.3-GHz  plane  wave  is  assumed,  with  the  same  set  of 
parameters  as  used  in  section  3.4.  The  medium,  9075  m  in  extent,  was  represented 
by  forty-nine  equally  spaced  phase  screens.  In  the  vertical,  a  sampling  interval  of 
0.2381  m  was  used  with  1024  points. 

The  FFT  routine  of  Press  et  al.^4  ^as  used  to  implement  the  split-step 
algorithm  of  Equation  4.1.  The  random  number  generator  RANI  developed  by 
Press  et  al.  was  used  in  the  random  phase  model  of  Appendix  A  to  produce  the 
phase  screens.  T1.&  autocorrelation  of  the  log-amplitude,  the  phase,  and  the 
scattered  field  of  each  realization  was  estimated  using  a  standard  unbiased 
estimator.35  The  procedure  was  repeated  for  250  independent  realizations:  the 


H.  Press,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W.  T.  VeUerling,  Numerical 
Recipes,  Cambridge  University  Press,  Cambridge,  England  (1986). 

V.  Oppenheim  and  R.  V.  Schafer,  Digital  Signal  Processing,  Prentice-Hall, 
Englewood  Cliffs,  N.  J.  (1975). 
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resulting  averaged  estimates  and  the  theoretical  prediction  are  plotted  in  Figure  4.1. 
The  autocorrelation  of  the  scattered  Held  displays  excellent  agreement  with  the 
theoretical  result  calculated  via  Equation  3.13.  At  zero  lag,  the  observed  variance 
of  the  log-amplitude  fluctuation  is  9.S6  x  10~^,  which  is  in  good  agreement  with 
the  predicted  result  (9.51  x  10"^)  in  section  3.4. 

As  expected,  the  two-dimensional  numerical  algorithm  accurately  predicts  the 
vertical  autocorrelation  of  the  three-dimensional  scattered  Held,  but  is  unable  to 
produce  the  desired  log-amplitude  fluctuations. 

HOMOGENEOUS  ISOTROPIC  TURBULENCE 

Within  the  constraints  discussed  in  section  4.1,  Equation  4.10  is  valid  for  any 
three-dimensional  homogeneous  spectrum.  The  numerical  example  in  section  3.2 
used  a  Gaussian  spectrum.  One  of  the  attractive  features  of  this  spectrum  is  that  it 
is  separable;  the  three-dimensional  spectrum  can  be  written  as  the  product  of  three 
one-dimensional  spectra.  This  simplifies  the  calculation  of  the  transverse 
spectrum.  Spectral  models  for  internal  waves  and  fme  structure  in  the  ocean  are 
also  separable.^^  Homogeneous  isotropic  turbulence  is  unfortunately  not 
separable,  and  the  resulting  transverse  spectrum  must  be  expressed  in  terms  of 
special  funedons.  In  this  secdon,  we  calculate  the  transverse  spectrum  and 
associated  autocorreladon  fimedon  for  homogeneous  isotropic  mrbulence. 

The  von  Karman  spectrum  is  used  to  describe  the  three-dimensional  index  of 
refracdon  fluctuadons.  Subsdtuting  Equadon  2.S  into  Equadon  4.10  yields 


St  ik, )  =  0.066(2rr)2  SxCl  axpi-kj  /  kI  ) 


dk^ 


JQ 


exp(-lj/Ar^) 

()k2+i2+^2)n/6  • 


(4.11) 


After  some  manipuladon,  the  integral  in  Equadon  4.1 1  can  be  expressed  exaedy  in 
terms  of  the  confluent  hypeigeometric  funedon 


(Kj  +  j 


■'(4.12) 


Equadon  4.12  can  be  simplified  in  the  inertial  subrange.  For  spadal  frequencies  k 
in  the  inerdal  subrange  and  below,  « 1 .  Hence  the  exponendal  can  be  set 
to  1  and  v^can  be  replaced  by  its  small  argument  approximadon, 


(4.13) 


E.  Ewart  and  S.  A.  Reynolds,  "Experimental  Ocean  Acoustic  Field  Moments 
versus  Predictions,"  in  Ocean  Variability  and  Acoustic  Propagation,  J.  Potter  and 
A.  Wam-Vamas  (ed.),  Kluwer  Academic  Publishers,  Boston,  Mass.,  pp.  2340 
(1991). 

^’M.  Abramowitz  and  I.  M.  Stegun,  Handbook  of  Mathematical  Functions,  Nauonal 
Bureau  of  Standards,  Washington,  D.  C.  (1965). 
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Autocorrelation 


Average  Autocorrelation  Estimate 

Gaussian  Spectrum 


Figure  4.1  Autocorrelations  of  field  that  has  propagated  through  a  medium  with  a 
Gaussian  spectrum.  Medium  of  variance  10-''2  and  correlation  length  3.04  m  are 
probed  by  3.3-GHz  plane  wave.  Theoretical  scattered  field  autocorrelation  is 
compared  to  numerical  result.  Phase  and  log-amplitude  autocorrelations  are  also 
shown.  Statistics  are  based  on  an  ensemble  average  over  250  realizations.  Each 
realization  had  forty-nine  equally  spaced  phase  screens  of  1024  points  and  a 
vertical  sampling  interval  of  0.2381  m.  Read  2.5e-4  as  2.5  x  1 0"^. 
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where  F  is  the  gamma  function.  For  parabolic  equation  simulations,  the 
maximum  spatial  firequency  of  interest  is  determined  by  the  wavelength  X  of  the 
illuminating  wave.  Since  typically  for  microwaves  X>tQ,  the  maximum 
frequency  divided  by  K„  wiU  be  much  less  than  1,  and  Equation  4.13  can  be  used 
for  all  spatial  frequencies  in  parabolic  simulations.  For  spatial  frequencies  in  the 
inertial  subrange,  we  observe  that  S,  is  proportional  to  , 

Corresponding  to  the  forward  transform  in  Equation  4.7,  the  transverse 
autocorrelation  function  is  given  by  the  inverse  transform 

=  dkMk2)txp{ikjZa).  (4.14) 

Substituting  the  approximate  transverse  spectrum  of  Equation  4.13,  the  integral 
can  be  expressed  in  terms  of  a  modified  Bessel  functiCHi,^* 

=  /io).  (4.15) 


Using  Equation  2.7,  the  transverse  autocorrelation  can  also  be  written  in  terms  of 
the  variance  of  the  index  of  refraction  <  >  as 


B,(rd)  =  2.647&(n}}io  '^^^576(^4  /^o).  (4.16) 


where  we  have  evaluated  the  gamma  function.  Equation  4.16  will  be  used  as 
"truth"  for  the  numerical  simulations.  This  factorization  is  convenient  for 
evaluating  by  a  series  expansion. 

We  note  that  Equation  4.16  is  indepen^nt  of  the  turbulent  inner  scale  fp.  For 
microwaves,  X  »  /q  and  the  probing  wave  does  not  "see"  the  very  small  scale 
fluctuations.  This  is  another  case  where  the  microwave  problem  differs  frotn 
analogous  propagation  problems  in  optics. 

AsZ(t->0,  the  small  argument  expansion  of  the  modified  Bessel  function  can 
be  used  to  yield  the  variance  of  the  transverse  process: 


B,iO)  =  lA94Sx{n})Lo, 


(4.17) 


The  transverse  variance  is  proportional  to  the  variance  of  the  medium  and  to  the 
outer  scale  of  the  turbulence  Lq. 

Individual  realizations  of  <j)t  were  generated  using  the  transverse  spectrum  of 
Equation  4.13  and  the  algorithm  in  Appendix  A.  Figure  4.2  shows  a  1024-point 


^*K.  H.  Craig  and  M.  F.  Levy,  "Parabolic  Equation  Modelling  of  the  Effects  of  Multipath  and 
Ducting  on  Radar  Systems,"  lEE  Proc.  F  138,  153-162  (1991). 
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realization  with  Lq  =  10  m,  step  size  Sx  =  100  m,  and  unit  variance  <nf>. 
Figure  4.3  shows  an  example  where  the  outer  scale  has  been  increased  to  100  m. 
This  example  vividly  exhibits  both  the  large-  and  small-scale  fluctuations 
characteristic  of  power-law  spectra.  For  comparison.  Figure  4.3  also  shows  a 
realization  generated  using  the  same  set  of  random  numbo's  and  the  same  outer 
scale  but  with  a  Gaussian  spectrum,  as  was  used  in  section  4.2.  In  contrast  to  a 
turbulent  process  that  has  both  an  inner  and  outer  scale,  a  Gaussian  spectrum  is 
characterized  by  a  single  correlation  length.  As  a  result,  the  realization  exhibits 
little  of  the  detailed  strucoire  present  in  turbulence  example. 

Five  hundred  realizations  of  the  process  with  Lo=  10  m  were  generated.  The 
autocorrelation  of  each  was  estimated  using  a  standard  unbiased  estimator  (see  Ref. 
35,  p.  539).  The  estimates  were  then  averaged  across  the  ensemble  and  compared 
to  the  theoretical  expression  in  Equation  4.21.  The  agreement  is  excellent  and  is 
shown  in  Figure  4.4.  The  study  was  repeated  with  the  outer  scale  increased  to  100 
m.  The  average  of  5000  estimates  are  compared  to  theory  in  Figure  4.5  with 
similar  results.  When  the  correlation  length  was  increased,  it  was  necessary  to 
increase  the  number  of  realizations  to  get  good  agreement 

Realizations  of  the  turbulent  process  with  Lq  =  10  m  were  used  in  a 
propagation  study.  Aside  from  the  different  spectrum  and  outer  scale,  the 
remaining  parameters  are  unchanged  from  the  numerical  example  in  section  4.2. 
Two  hundred  fifty  independent  simulations  were  used  to  estimate  the  ensemble 
statistics  of  the  field.  Figure  4.6  compares  the  predicted  autocorrelation  of  the 
scattered  field  to  three-dimensional  theory.  The  agreement  is  excellent,  particularly 
at  short  lags.  Beyond  the  correlation  length  of  the  medium,  the  error  becomes 
more  noticeable.  This  is  to  be  expected,  since  less  information  is  available  at  large 
lags  to  make  the  estimates.  To  improve  the  agreement,  the  averaging  could  be 
conducted  over  more  realizations.  But  since  the  validity  of  the  turbulence  model 
itself  becomes  questionable  beyond  the  correlation  length  (section  2),  it  is  probably 
not  worth  the  effort. 

Also  shown  in  Figure  4.6  are  the  estimated  log-amplitude  and  phase 
autocorrelation.  Since  the  two  curves  differ  significantly,  we  are  clearly  not  in  the 
two-dimensional  Fraunhofer  regime.  At  zero  lag,  the  simulations  give  the  log- 
amplitude  variance  to  be  1.02  x  10“^.  The  theoretical  value  calculated  by  Equation 
3.9  is  3.19  X  1(H.  This  example  again  illustrates  that  while  the  cross-sectional 
model  produces  good  results  for  the  moments  of  the  field,  it  fails  to  calculate  the 
log-ampliuide  fluctuations. 

4.4  SURFACE  LAYER  TURBULENCE 

The  procedure  for  generating  random  realizations  must  be  modified  in  the 
aurospheric  surface  layer.  Sections  2.3  and  2.4  showed  that  at  low  altitudes  both 
the  sttucuire  constant  cl  and  tlie  outer  scale  Lq  must  be  treated  as  height-dependent 
variables.  The  only  study  known  to  model  electromagnetic  propagation  through 
surface  layer  turbulence  using  the  parabolic  equation  was  conducted  by  Levy  and 
Craig.^'^®  In  this  section,  the  algorithm  used  by  Levy  and  Craig  is  critiqued,  and 
it  is  suggested  that  they  used  an  incorrect  power  law  specttum  in  their  simulations. 
Possible  modifications  to  existing  meteorological  models  for  surface  layer 
turbulence  are  then  considered  to  produce  expressions  compatible  with  our  simula- 
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FIgurs  4.3  Realization  of  Same  as  Figure  4.2  except  with  outer  scale 
increased  to  100  m.  For  comparison,  realization  of  a  process  with  the  same  outer 
scale  but  with  a  Gaussian  spectrum  is  also  shown.  Step  size  is  equal  to  100  m. 
Unit  variance  and  sampling. 
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Figure  4.4  Transverse  autocorrelation  for  process  with  Lo>*  '[O  m.  Theory  (solid 
curve)  compared  to  ensemble  average  of  500  realizations  (dots),  using  1024  point 
realizations.  Step  size  is  equal  to  100  m. 
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Average  Autocorrelation  Estimate 

Kolmogorov  Spectrum 


Figure  4.6  Autocorrelations  of  field  that  has  propagated  through  a  medium  with  a 
Kolmogorov  spectrum.  Medium  of  variance  10”’^  and  outer  scale  10  m  are  probed 
by  3.3-GHz  plane  wave.  Theoretical  scattered  field  autocorrelation  is  compared  to 
nu.merical  result.  Phase  and  log-amplitude  autocorrelations  are  also  shown. 
Statistics  are  based  on  an  ensemble  average  over  250  realizations.  Each 
realization  had  forty-nine  equally  spaced  phase  screens  of  1024  points  and  a 
vertical  sampling  interval  of  0.2381  m.  Read  7.0e-4  as  7  x  lO"^. 


41 


tion  approach.  Finally,  the  modified  spectra  are  used  to  calculate  the  transverse 
spectra  necessary  for  parabolic  equation  simulations. 

4.4.1  Approach  of  Levy  and  Craig 

The  randcnn  index  of  refraction  nyis  inhomogeneous  and  possibly  anisotropic  in 
the  atmospheric  surface  layer.  Unfortunately,  the  spectral  filtering  method  for 
generating  realizations  is  limited  to  homogeneous  processes  and  hence 
inapplicable.  To  circumvent  these  difficulties.  Levy  and  Craig  define  an 
intermediate  process  q  =  nf  I C^.  They  then  make  the  simplifying  assumption 
that  ^  is  statistically  homogeneous  and  isotropic  and  therefore  can  be  simulated 
using  spectral  filtering.  By  assuming  Sx  »Lo,  they  generate  independent  one¬ 
dimensional  realizations  of  ^  at  uniform  range  steps.  They  use  the  one¬ 
dimensional  form  of  the  Kolmogorov  spectrum  that  shows  a  fiequency 
dependence  to  design  the  spectral  filter  (Eq.  2.1)  and  then  the  produce  realizations 
using  the  algorithm  in  Appendix  A.  The  fluctuating  index  of  refraction  is 
recovered  by  rt/  =  C„|.  Numerical  propagation  studies  are  presented  without 
comparison  to  theory. 

The  main  error  made  by  Levy  and  Craig  is  in  selecting  the  wrong  spectrum  to 
generate  realizations.  The  quantity  of  relevance  to  parabolic  simulations  is  the 
integrated  refractive  index  fluctuations  (Eq.  3.4b);  consequently,  realizations  should 
be  generated  using  the  transverse  spectrum  S,  (Eq.  4.7).  For  turbulent  fluctuations. 
Si  exhibits  a  frequency  dependence  in  the  inertial  subrange  rather  than  the 
power  law  used  by  Levy  and  Craig.  Selecting  the  wrong  spectrum  also 
introduces  errors  in  scaling  and  dimensionality.  (The  one-dimensional  spectrum 
has  dimensions  of  meters,  whereas  St,  like  S„,  has  dimensions  of  meters  cubed.) 

4.4.2  Modified  Spectral  Models  for  Surface  Layer  Turbulence 

Levy  and  Craig  define  an  (incorrect)  intermediate  process  and  then  assume  it  is 
stationary.  A  more  fundamental  approach  is  to  return  to  the  three-dimensional 
spectral  model  for  the  turbulence  and  then  determine  what  properties  it  must 
possess  to  force  an  appropriately  defined  intermediate  process  to  be  stationary.  In 
this  section,  two  possible  models  for  boundary  layer  turbulence  are  considered  that 
might  be  numerically  tractable.  Both  models  treat  the  specDul  bandwidth  as  a 
height-independent  constant  One  apjHoach  considers  the  height-dependent  strucuire 
constant  to  be  a  separable,  deterministic  function.  The  second  model  treats  the 
turbulence  as  a  quasi-nomogcncous  pro^-ss.  In  section  4.4.3,  the  transverse 
spectra  for  both  models  are  discussed. 

We  begin  with  the  standard  one-dimensional  spectrum  for  inhomogeneous, 
height-dependent,  surface  layer  turbulence  ^ 


V,{K-,z)  =  0.mci{z)K-^<\  (4.18) 

where  the  strucoire  constant  can  be  written  in  terms  of  the  turbulent  outer  scale 
and  the  gradient  of  the  mean  profile  (Eq.  2.9).  For  wave  propagation  studies,  the 
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full  three-dimensional  spectrum  is  required.  Assuming  the  medium  is  isotropic, 
the  three-dimensional  Kolmogorov  specuum  follows  from  Equation  2.1; 

5„(A::z)  =  0.033C2(2)^:-"/1  (4.19) 


Although  it  is  not  often  mentioned  in  the  atmospheric  literature,  spectral 
representations  such  as  those  given  in  Equation  4.19  are  limited  to  frequencies  in 
the  inertial  subrange.  From  section  2,  the  inertial  subrange  is  defined  by  spatial 
frequencies  within  the  band  lit!  Lq  <K  <lnl  Iq.  In  the  surface  layer,  however, 
the  outer  scale  Lq  is  itself  a  function  of  height  through  Equation  2.10.  Hence  we 
have  the  cumbersome  situation  where  both  the  amplitude  (through  C^)  and  the 
lower  break  frequency  {lit! Lq)  of  the  surface  layer  spectrum  are  dependent  on 
height.  For  example,  inferring  an  inhomogeneous  von  Karman  specu'um  directly 
from  Equations  2.5  and  4.19  yields 


5„(K:2)  =  0.033C2(2)  exp(-K2/A:2)[(^2^^2(^)J-ll/6 

Following  the  procedure  of  section  4.3,  an  inhomogeneous  uansverse  spectrum  can 
presumably  be  derived.  But  since  the  lower  break  frequency  of  the  spectrum  will 
depend  on  height,  realizations  of  the  uansverse  process  cannot  be  generated  using 
the  filtered  white  noise  approach  detailed  in  Appendix  A. 

As  a  first  approximation,  we  retain  the  height  dependence  through  but 
replace  the  explicit  Lq{z)  in  Equation  4.20  with  a  height-independent  constant  Lc. 
This  is  probably  reasonable,  since  we  arc  primarily  interested  in  the  inertial 
subrange  of  turbulence  and  in  this  regime  the  only  height  dependence  appears  in 
Cn-  In  making  this  approximation,  the  spectral  bandwidth  of  the  inertial  subrange 
is  made  height  independent.  To  estimate  an  appropriate  numerical  value  for  L^ 
one  possible  approach  is  to  relate  it  to  the  variance  of  the  index  of  refraction.  This 
is  detailed  in  Appendix  B. 

There  is  a  second  problem  with  the  hypothesized  spectrum  in  Equation  4.20. 
This  inhomogeneous  specpum  can  presumably  be  related  to  a  sUucture  function  or 
an  autocorrelation  by  an  appropriate  uansform  pau.  These  are  two-point  statistical 
quantities,  but  Equation  4.20  depends  on  only  a  single  height  z.  A  more  general 
model  would  show  the  explicit  two-point  dependence. 

With  these  two  conditions  in  mind,  we  hypothesize  a  modified  three- 
dimensional  spectral  model  for  turbulence  as 


(4.21) 


where 


0„(K)  =  cxp(-/('2/k2  )(/!:2  4-  Z.-^r>  ,  (4.22) 
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This  model  partitions  the  turbulence  into  a  height-independent  spectrum  <1>„  and 
the  quantity  C  (21.22).  which  in  general  depends  on  two  heights.  There  is  some 
latitude  in  defining  but  it  should  be  chosen  to  be  consistent  with  existing 
meteorological  expressions.  We  will  consider  two  possible  models.  In  the  first, 
the  structure  constant  is  treated  as  a  deterministic  quantity.  The  ambiguous  term 
(z)  is  separated  to  show  explicit  two-point  dependence: 


C(zi.Z2)  =  0.033C„(z,)C„(z2).  (4.23a) 


Note  that  at  zi  =  Z2,  Equation  4.23a  is  consistent  with  the  standard  representations 
in  Equations  4.18  and  4.19. 

A  second  model  treats  the  medium  as  being  quasi-homogeneous.  The 
two-point  statistics  of  a  quasi-homogeneous  process  can  be  written  as  the  product 
of  a  term  depending  on  the  difference  coordinate  and  a  term  depending  on  the 
average  position.  The  spectrum  <I>„  is  the  Fourier  transform  of  the  difference 
coordinate  component,  and  hence  C  is  written  as 


C(z„22)  =  0.033C2[(z,  +Z2)/2].  (4.23b) 


The  component  of  a  quasi-homogeneous  process  that  depends  on  the  average 
coordinate  is  generally  a  slowly  varying  function  of  its  argument.  Turbulent  media 
are  often  modeled  as  quasi-homogeneous  processes  in  theoretical  wave  propagation 
studies.'^'’'^® 

For  either  model,  it  is  useful  to  define  the  autocorrelation  of  the  medium  by  an 
inverse  transform. 


(nyr(xi,yi,z,)ny:(x2,y2.22))  = 

-C(2i.22)r  dK«I>(K)exp(iKT),  (4.24) 

J  —OP 

where  K  =  {k^,  ky,  k^)  and  r  =  {xd.yd.^d)-  This  formulation  is  u.seful  for 
calculating  the  transverse  spectrum  for  inhomogeneous  media.  The  .'.iverse 
transfomi  is  evaluated  in  Appendix  B. 

4.4.3  Transverse  Spectrum 

The  pertinent  statistical  quantities  for  parabolic  equation  simulations  are  the 
autocorrelation  of  the  integrated  refractive  index  fluctuations  0,  and  the  associated 
transverse  specUum.  Using  Equations  4.4b  and  4.24,  and  again  assuming  d-.e  step 
size  is  large  compared  to  the  correlation  length,  it  follows  that 
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=  &cC(zi,Z2)J  dl^  j  dk'^ik^  =  0,ky,lc^exp(ikjZj) . 

J  — «o  J  — «»  / 


(4.25) 


Equation  4.25  is  valid  for  either  form  of  We  now  consider  the  separable 

model  of  Equation  4.23a.  It  is  useful  to  define  the  intermediate  function, 

^(2)  =  (O.O33)-*/20,(2)/C„(z).  (4.26) 


Since  C  is  deterministic  in  this  model,  it  can  be  brought  inside  the  ensemble 
averaging  brackets  and 


($(2i)^(22)) = &  r  dky  r 

J  «•»  J  —00 


dk^<t>(kjc=0,ky,k^  txpi+ikizj) 


mB^(zi-Z2). 


(4.27) 


The  autocorrelation  of  |  depends  solely  on  die  difference  coordinate;  hence  the 
intermediate  process  is  statistically  homogeneous.  Taking  the  Fourier  transform 
gives  the  spectrum  S^: 


r(ii/6) 


(4.28) 


This  spectrum  can  be  used  with  the  algorithm  in  Appendix  A  to  generate 
realizations  of  Realizations  of  <l>(  that  are  necessary  for  the  parabolic  equation 
simulations  are  recovered  via  Equatiai  4.26. 
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5.  SUMMARY  AND  SUGGESTIONS  FOR 
FURTHER  WORK 


Based  on  the  analysis  and  numerical  results  of  the  preceding  sections,  we  can 
conclude  that  the  parabolic  equation/split-step  algorithm  is  a  viable  approach  for 
modeling  electromagnetic  propagation  through  atmospheric  turbulence.  The 
propagation  studies  presented  in  this  report  were  generated  independently  from 
TEMPER.  It  is  a  relatively  simple  task  to  design  a  subroutine  that  would  also  add 
"randomness"  to  the  existing  TEMPER  code;  all  that  is  required  is  an  appropriate 
spectral  model,  a  random  number  generator,  and  access  to  a  fast  Fourier  transform 
routine.  The  procedure  should  produce  reliable  results  in  an  ensemble  sense;  the 
correct  second-order  statistics  of  the  propagating  field  can  be  predicted  by  averaging 
across  many  realizations.  As  used  in  this  report,  spectral  modeling  assumes 
Gaussian  statistics  for  the  random  medium.  Consequently,  the  method  is  incapable 
of  modeling  low-probability  events  governed  by  non-Gaussian  statistics. 

Tlie  specific  conclusions  are  summarized  as  follows: 

1.  Although  the  parabolic  equation  method  has  previously  been  used  to  study 
optical  and  acoustical  propagation  through  random  media,  microwave  propagation 
presents  a  unique  set  of  problems.  Since  the  perturbations  to  the  field  are 
relatively  small,  weak  scattering  theory  can  usually  be  used  to  estimate  the 
propagating  field.  The  inhomogeneous  and  possibly  anisouopic  nature  of  surface 
layer  turbulence  must  be  considered  for  low-altitude  propagation. 

2.  Two-dimensional  models  like  TEMPER  are  probably  adequate  for 
calculating  the  statistical  moments  of  the  propagating  field.  They  may  not, 
however,  suffice  for  calculation  of  the  log-amplitude  fluctuations  except  in  the  very 
far  field.  The  term  "far  field"  itself  has  different  interpretations  in  two  and  three 
dimensions. 

3.  Provided  that  the  step  size  is  targe  compared  to  the  turbulent  outer  scale,  the 
medium  can  be  modeled  by  uncorrelated  random  phase  screens.  The  relevant 
statistical  quantity  for  generating  the  phase  screens  is  the  transverse  spectrum.  For 
homogeneous  isotropic  turbulence  in  the  inertial  subrange,  the  transverse  spectrum 
shows  a  frequency  dependence.  This  is  in  conuast  to  the  one-dimensional 
spectrum  with  a  power  law  used  by  Levy  and  Craig. 

4.  The  transverse  spectrum  of  homogeneous  isotropic  turbulence  can  be 
expressed  in  terms  of  a  confluent  hypergeometric  function  (Eq.  4.12).  At 
microwave  frequencies,  a  small  argument  approximation  can  be  used  to  show  the 
/f-8/3  behavior  in  the  inertial  subrange  (Eq.  4.13).  For  a  statistically 
homogeneous  medium,  the  transverse  autocorrelation  can  be  evaluated  and 
expressed  using  a  modified  Bessel  function  (Eq.  4.15). 

5.  Turbulence  is  inhomogeneous  in  the  surface  layer  comprised  of  altitudes  less 
than  approximately  100  m.  Both  the  structure  constant  and  the  outer  scale  depend 
on  height.  Since  the  outer  scale  determines  the  lower  spatial  frequency  bound  of 
the  inertial  subrange  (Fig.  3.1),  the  implication  is  a  specu-um  where  both  the 
amplitude  and  the  bandwidth  depend  on  altitude.  Such  a  spectrum  is  incompatible 
with  the  parabolic  equation/split-step  algorithm.  As  an  approximation,  the 
bandwidth  was  made  constant  and  independent  of  height.  Two  plausible  forms  of 
the  structure  constant  were  considered.  In  one  model,  the  structure  constant  was 


46 


treated  as  a  separable,  deterministic  function.  In  the  second  model,  the  structure 
constant  was  treated  as  a  quasi-homogeneous  function. 

The  expressions  derived  for  the  transverse  spectra  and  autocorrelations  derived  in 
section  4  are  amenable  to  numerical  implementation.  A  limited  number  of 
independent  simulations  were  presented  in  this  report.  A  possible  program 
validating  the  proposed  algorithm  when  coupled  with  TEMPER  is  outlined  below: 

1.  Classical  results  for  the  canonical  problems  should  be  used  as  benchmarks 
in  validating  the  computer  routines.  Analytical  solutions  are  available  for  point 
source  and  plane  wave  propagation.  Since  these  pristine  solutions  neglect 
interaction  with  deterministic  refractive  index  profiles  and  rough  surfaces,  these 
features  should  be  "turned  off  when  first  testing  randomness  in  TEMPER. 

2.  Initial  numerical  tests  for  both  the  plane  wave  and  point  source  problems 
might  best  be  done  by  assuming  a  Gaussian  spectrum  for  the  medium.  Because  a 
Gaussian  spectrum  is  characterized  by  a  single  scale  size,  the  resulting  realizations 
do  not  exhibit  the  detailed  fine  structure  of  processes  with  power-law  spectra  (Fig. 
4.3).  Consequently,  Monte  Carlo  simulations  should  converge  relatively  rapidly.. 
Agreement  with  the  simulations  presented  in  section  4.2  should  be  expected. 

3.  The  second  set  of  tests  for  both  the  plane  wave  and  the  point  source  can  be 
conducted  for  homogeneous  isotropic  turbulence.  The  effect  of  varying  the 
numerical  values  of  the  various  parameters  (outer  scale,  "universal  constants,"  etc.) 
can  be  studied.  Required  sampling  rates  can  be  defined,  and  the  number  of 
realizations  necessary  for  convergence  of  the  ensemble  averages  can  be  estimated. 
The  numerical  experiments  pre.sented  in  section  4.3  should  be  replicated. 

4.  Finally,  the  full  features  of  TEMPER  can  again  be  "turned  on."  The 
structure  constant  for  representative  atmospheric  profiles  can  be  calculated. 
Realizations  of  inhomogeneous  surface  layer  turbulence  can  be  generated  and 
included  in  TEMPER.  The  variance  of  the  propagating  field  can  be  estimated. 
Numerical  results  should  ultimately  be  compared  to  experimental  observations. 

Beyond  the  simple  expedient  of  adding  atmospheric  randomness  to  TEMPER, 
this  work  has  suggested  directions  for  continued  research: 

1.  Independent  from  TEMPER,  additional  simulations  should  be  done  for  the 
point  source  problem.  The  analysis  comparing  the  two-  and  three-dimensional 
problems  should  also  be  repeated  for  this  configuration. 

2.  A  good  three-dimensional  spectral  model  for  boundary  layer  turbulence  is 
needed.  For  a  lack  of  quantitative  data,  the  emphasis  in  this  report  is  on  isotropic 
models.  Boundary  layer  turbulence  is  surely  anisotropic,  with  relatively  long 
correlation  lengths  in  the  horizontal.  It  would  also  be  of  interest  to  analyze  the 
effect  anisotropy  would  have  on  the  various  statistics  of  the  field. 

3.  Numerical  techniques  could  be  developed  for  simulating  quasi-homogeneous 
random  media.  It  seems  likely  that  an  algorithm  could  be  developed  that  would 
parallel  that  used  for  homogeneous  media  in  Appendix  A.  Such  an  algorithm 
would  also  be  of  interest  to  researchers  working  in  rough  surface  scattering. 

4.  Fractals  appear  to  be  a  promising  method  for  simulating  turbulence;  their 
most  attractive  feature  is  the  ability  to  model  intcrmiiicncy.  A  good  appreciation 
of  atmospheric  dynamics  is  needed  to  understand  the  role  of  inicrmiuency  and  how 
it  should  be  modeled  with  fractals.  Further  analytical  work  is  necessary  to  predict 
the  properties  of  a  wave  propagating  through  a  fractal  medium. 

5.  To  validate  an  upgraded  version  of  TEMPER  containing  turbulent 
fluctuations,  the  philosophy  of  this  report  has  been  to  compare  the  numerical 


results  to  established  theoretical  solutions.  Only  theoretical  solutions  to  relatively 
simple  problems  have  been  considered,  such  as  plane  wave  and  point  source 
propagation.  Solutions  to  other  problems  are  available  in  the  literature;  for 
example,  Ishimaru^  considers  the  propagation  of  a  Gaussian  beam.  Path  integral 
methods  that  reportedly  combine  large-scale  deterministic  features  with  randomness 
might  be  adapted  to  our  purposes.^^-^®  Modal  techniques  for  propagation  through 
deterministic  profiles  perhaps  could  be  modified  to  include  randomness. 

6.  At  microwave  frequencies,  the  turbulence  possesses  structure  on  the  same 
scale  as  the  wavelength.  It  is  possible  that  this  would  induce  polarization  effects 
that  are  ignored  in  our  scalar  model.  If  possible,  the  magnitude  of  these  effects 
should  at  least  be  estimated.  It  is  unclear  whether  these  features  could  be  included 
in  TEMPER. 
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APPENDIX  A. 


GENERATING  RANDOM  REALIZATIONS 

Realizations  of  atmospheric  turbulence  can  be  simulated  in  various  ways.  The 
most  complete  method  is  to  solve  the  governing  Navier-Stokes  equations  describing 
the  fluid  motion  numerically.  Unfortunately,  this  procedure  is  not  feasible  at  the 
large  scales  necessary  for  wave  propagation  studies.  As  an  alternative,  we  are  left 
with  numerical  techniques  that  produce  realizations  consistent  with  some  order 
(usually  second)  of  the  statistics  describing  the  medium.  While  these  methods  yield 
realizations  exhibiting  the  correct  statistics,  the  individual  realizations  may  or  may 
not  actually  "look"  like  the  true  process.  As  a  result,  these  methods  are  best  used  as 
part  of  an  algorithm  where  the  entire  procedure  (in  this  case,  wave  propagation 
through  a  random  medium)  is  simulated  many  times  and  then  the  statistics  are 
inferred  by  averaging  across  an  ensemble  of  realizations.  TTiese  procedures  fall  under 
the  general  heading  of  Monte  Carlo  techniques. 

One  of  the  first  studies  of  wave  propagation  through  random  media  using  the 
parabolic  equation  was  conducted  by  Macaskill  and  Ewart. ''*  Following  their 
example,  most  subsequent  studies  have  simulated  realizations  of  the  random  medium 
using  spectral  filtering.  In  this  approach,  realizations  of  Gaussian  white  noise  are 
first  generated  using  a  random  number  generator  and  then  filtered.  To  obtain  the 
desired  output  process,  the  filter  is  designed  based  on  the  square  root  of  the  desired 
spectrum.''^  An  efficient  algorithm  based  on  the  fast  Fourier  transform  (FFT)  is  now 
presented.  This  approach  is  particularly  attractive  in  the  present  context,  since  the 
FFT  is  already  used  in  the  parabolic  equation/split-step  algorithm. 

Consider  an  N  point  process  0,(^J  defined  by 

t  =  0X...N-h  (Al) 


where  the  square  brackets  indicate  a  sampled  version  of  the  continuous  process  and 
the  sampling  interval  is  8z.  The  discrete  spectrum  is  a  sampled  version  of  the 
continuous  spectrum 


S, [;■]  =  (<52)-' S, {j/{N5z)l  j  =  0X...  N 12 

,  (A2) 

S,  [7]  =  (<5z)-'  S,  {{N-j)/{Ndz)),  j  =  N/2, 

The  discrete  realization  is  recovered  via  an  inverse  discrete  transform, 

N-l  _ 

=  S  0,(7)  exp[-/27ry//iV],  (A3) 

j=0 

where  the  spectral  components  are  generated  according  to  the  rule 

0([7]  =  [-'S([7]A' ln(<7y)]'^^exp(t2;rrp  for  jtO,j^NI2 

(A4) 

0([7]  =  [-'^(t7]2A^  ln((?p]'^^cos(2w^)  for  7  =  0  or  j  =  NI2. 


*’C  Macaskill  and  F  1:  Kwart,  "Compuier  Simulation  of  I  wo  Dimensional  Random  Wave  Propagation," 
IMA  J  Appl  Math  33.  1-15  (1984) 

Papoulis,  Probability,  Random  VartaoUs,  and  Stochastic  Processes,  McGraw-Hill,  New  York  (1984) 
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In  Equation  A4,  and  are  random  numbers  uniformly  .istributed  between  zero  and 
1.  Note  that  both  the  zero  frequency  and  the  Nyquist  component  are  purely  real.  To 
make  the  realization  purely  real,  there  is  the  additional  symmetry  constraint  that 


^ “  ^N—J  *  j  t ,  N I  'I 

rj=rf^-j,  j  =  .N 12. 


(A5) 


The  inverse  transform  in  Equation  A3  is  efficiently  evaluated  using  an  FFT.,  The 
resulting  realization  of  the  process  is  purely  real.  A  marginal  increase  in  efficiency 
can  be  gained  by  relaxing  the  symmetry  constraint,  generating  a  complex  <p,,  and 
then  treating  the  real  and  imaginary  parts  as  independent  realizations.^^  Explicit 
evaluation  of  the  probability  integrals  shows  that  this  algorithm  produces 
realizations  with  the  desired  spectrum. 


M  Minin  and  S.  M.  Flute,  ■Intensity  Images  and  Siaiistics  from  Numerical  Simulation  of  Wave 
Propagauon  in  3-D  Random  Media,"  Appl  Opt  27,  2111-2127  (1988) 
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APPENDIX  B. 


VARIANCE  OF  REFRACTIVE  INDEX  IN  THE 
SURFACE  LAYER 

As  discussed  in  section  atmospheric  turbulence  is  typically  modeled  as  a 
locally  homogeneous  random  process;  as  such,  it  must  be  characterized  by  a  structure 
function  rather  than  an  autocorrelation.  In  the  frequency  domain,  the  locally 
homogeneous  character  is  reflected  in  the  spectrum  being  undefined  for  the  low 
spatial  frequencies  in  the  energy  input  regime.  For  mathematical  convenience, 
however,  a  spectral  shape  is  often  assigned,  making  an  inverse  transform  possible, 
and  hence  a  function  that  can  be  interpreted  as  an  autocorrelation. 

In  section  4.4.3,  a  model  for  surface  layer  turbulence  was  developed.  A  parameter 
Lj  was  defined  with  units  of  meters,  where  the  upper  bound  of  the  energy  input 
regime  was  defined  by  2nlLc.  One  method  for  assigning  a  numerical  value  to  is  to 
examine  its  relationship  to  the  variance  of  the  medium.  To  calculate  the  variance, 
we  must  first  calculate  the  autocorrelation. 

The  autocorrelation  was  related  to  the  spectrum  in  Equation  4.24.  If  the  spectrum 
is  isotropic,  the  three-dimensional  integral  can  be  converted  to  a  single  integral,®* 
and 


('4 :  zi  •  ^2 )  =  ?(zi .  Z2  )4  ro-J"* 


n 

Jo 


<i>{K)K&m{Kri)dK. 


(Bl) 


Let  2S2^=^2.  Substituting  the  spectrum  into  (Bl),  neglecting  the  exponential 
term,  and  evaluating  the  resulting  integral  yields®^ 


B,(rrf;z)  =  0.391C^(r)Lf^rrfL;’/2)''^Ari/3(rrfL-').  (B2) 


The  variance  is  derived  by  evaluating  Equation  B3  at  rj  =  0: 


=  B,(r<,=0)=  0.523C,^z)Z,p .  (B3) 


Hence  the  variance  of  the  medium  is  proportional  to  The  structure  constant 

C^(z)  is  given  by  Equation  2.9. 


^*A  Ishimani,  Wave  Propagation  arxt  Scattering  in  Random  Media,  Academic  Press,  New  York,  p  517 
(1978) 

S  Gradshicyn  and  I  M  Ryzhick,  Tables  of  Integrals,  Series  and  Products,  Academic  Press,  Orlando, 
Fla  (1980) 
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